Semiclassical approximation for the second harmonic generation in nanoparticles
Abstract
Second harmonic generation by spherical nanoparticles is a nonlocal optical process that can also be viewed as the result of the nonlinear response of the thin interface layer. The classical electrodynamic description, based e.g. on the nonlinear Mie theory, entails the knowledge of the dielectric function and the surface nonlinear optical susceptibility, both quantities are usually assumed to be predetermined, for instance from experiment. We propose here an approach based on the semiclassical approximation for the quantum sumoverstates expression that allows to capture the secondorder optical process from first principles. A key input is the electronic density, which can be obtained from effective single particle approaches such as the densityfunctional theory in the local density implementation. We show that the resulting integral equations can be solved very efficiently rendering thus the treatment of macroscopic systems. As an illustration we present numerical results for the magic Na cluster.
pacs:
78.67.n, 73.20.Mf, 71.10.w, 42.65.Ky, 31.15.A,73.22.DjI Introduction
For the description of a wide range of phenomena in physics and chemistry one is faced with the question of how to predict in a reliable and system specific way the response to external electromagnetic fields that impart energy and possibly momentum to the system Mahan (2000). To name but few examples, we mention here the response of nanoparticles to light in the linear Halas et al. (2011) and the nonlinear regimes Roke (2009). For electrons as a perturbation we refer to the overview García de Abajo (2010). A wellstudied route to address these issues is the concept of the dielectric response functions which is detailed in standard books such as in Mahan (2000), but also in recent reviews, e.g. in Refs. Halas et al. (2011); Roke (2009); García de Abajo (2010).
For extended bulk metals and metallic surfaces the dielectric response was extensively studied (Ref. Liebsch (1997); Andersen and Hübner (2002) and references therein). An extension to finitesize systems brings about a number of new aspects that complicate the theory but at the same time offer new opportunities for the occurrence of new phenomena. For example, it is known that the second harmonic generation (SHG) is forbidden in systems with an inversion symmetry on the dipole approximation level. Therefore Östling, Stampfli, and Bennemann Östling et al. (1993) and Dewitz, Hübner, and Bennemann Dewitz et al. (1996) proposed the anharmonic oscillator model to describe the secondorder nonlinear response from spherical systems. A further step was undertaken by Dadap et al. Dadap et al. (1999) who considered the smallparticle limit () and described SHG as a mixture of dipole and quadrupole excitation processes ( stands for the system size). It was also shown how to connect the model output to experimentally observed quantities, i.e. the elements of the nonlinear optical susceptibility tensor . In this approach the inversion symmetry is not a necessary assumption as it is assumed that the second harmonic generation originates from the thin surface layer, where the symmetry is broken anyway. This idea can be extended in several directions as was demonstrated in numerous works: Pavlyukh and Hübner Pavlyukh and Hübner (2004) developed the nonlinear Mie theory valid for particles of arbitrary sizes, de Beer and Roke included the sumfrequency generation mechanism into the considerations de Beer and Roke (2009), the cylindrical geometry was treated by Dadap Dadap (2008) and finally the theory for arbitrarily shaped particles was developed by de Beer, Roke, and Dadap de Beer et al. (2011).
All these theories based on classical electrodynamics rely on the knowledge of the frequencydependent dielectric function and the nonlinear optical susceptibility tensor. With the fabrication processes being perfected and the system sizes tending smaller and smaller one may wonder to which extent quantum effects are important and whether it is justified to use the same susceptibility tensor to describe semiinfinite and finite size systems on the nanometer range. To shed light on these issues, it is desirable to have a quantum theory for the nonlinear response on the nanoscale. The fully atomistic approach seems to be out of reach for present computers, as currently maximally hundreds of atoms are possible to treat using quantum chemistry codes. Yet, the outstanding question is, how important are the electronic correlation effects and is there possibly a way to stay on the solid quantum theory basis while treating larger systems?
There is an affirmative answer to these questions as was demonstrated recently in the linear optics case by Prodan and Nordlander Prodan and Nordlander (2003). They succeeded to push the limits of the timedependent density functional theory (TDDFT) to metallic systems containing millions of atoms. But at the same time they demonstrated that for these system sizes the semiclassical approach becomes very accurate. This is a marked finding as it allows to replace the complicated sumoverstates quantum mechanical expression Orr and Ward (1971) for the polarizability tensor with a single integral equation. Consequently, there is only one parameter in the theory: the ionic density distribution. With some reasonable assumptions about the ionic density such as in the jellium model (this assumption is reasonable even for molecular structures, as we have shown recently Pavlyukh and Berakdar (2009, 2010) for fullerenes. The usefulness of the jellium model was demonstrated by the pioneering works of Ekardt on sodium clusters Ekardt (1984, 1985) or of Puska and Nieminen on C Puska and Nieminen (1993).) we can obtain the ground state electronic density from the solution of the KohnSham equations and express the response function in its terms. The Drude dielectric function follows automatically.
The semiclassical approximation is rooted in the work of Mukhopadhyay and Lundqvist Mukhopadhyay and Lundqvist (1975) who derived the corresponding integral equation within linear response theory. Their theory was applied in numerous cases ranging from plasmons in the jellium model to collective resonances in C Vasvári (1996) or carbon nanotubes Vasvári (1997). The equation can also be derived starting from the quantum mechanical sumoverstates expressions Pavlyukh et al. (2012) and using the assumption that the frequency of the external field is large compared with the singleparticle gap () Ichikawa (2011).
One may wonder why this program was not implemented for nonlinear optical processes. As a matter of the fact, already in 1973 Wang, Chen and Bower Wang et al. (1973) classically treated second harmonic generation from alkali metals. A decade later Apell Apell (1983) derived an expression for the second harmonic current in the form:
(1) 
where for the unperturbed ground state electron density in the form the parameter is a function of and , whereas and are constants. The theory gained less attention than, for example, the Mukhopadhyay and Lundqvist work because the connection to quantum mechanics was lost. Here we mention nonetheless numerous works in the field extending over more than three decades: the small homogeneous spherical particle limit Hua and Gersten (1986), the RayleighGans scattering approximation for a lattice of such particles Martorell et al. (1997), second harmonic generation by twodimensional particles Valencia et al. (2003), or a more recent treatment of arbitrary geometries Zeng et al. (2009). Notice that the assumption of homogeneous electron density distribution within the sample is an inevitable component of such classical theories.
It is possible to revive the theory by noting that the electric field () in Eq. (1) must also include the induced field. This simple observation immediately raises the level of the theory to the random phase approximation (RPA) or, if we include electronic correlations, to TDDFT level. Our manuscript is, hence, a formalization of this message.
To this end we first derive an expression analogous to (1) starting from the sumoverstate quantum mechanical formula for the nonlinear optical susceptibility Orr and Ward (1971) and employing the highfrequency approximation. Although quite technical, we believe that this derivation (Appendix A) has its own merits as it establishes the equivalence of the hydrodynamic approximation of Apell and the highfrequency semiclassical expansion. It is interesting to notice that the asymptotic behavior of this quantity is at variance with the results of Scandolo and Bassani Scandolo and Bassani (1995) who predict a decay. This seems, however, to be a direct consequence of the assumption of the alldipole excitation mechanism that underly their study. It is possible to prove from general principles that the inclusion of the quadrupole excitations leads to the asymptotics Satitkovitchai (2009).
Our derivation raises the question of whether it is sufficient to know the unperturbed ground state density to obtain the lowest order approximation for an arbitrary response function. We recall that from the point of view of the diagrammatic perturbation theory Ward (1965) SHG comprises three processes in which the photon is emitted before, between or after the absorption of two photons. Consequently, one might wonder if each diagram of this expansion can also be expressed in terms of . As we demonstrate below, the answer is negative, one additionally needs the oneparticle density matrix whose diagonal elements are given by . This comes not as a surprise if we consider the analogy with the orbitalfree kinetic density functional theory Wang and Carter (2002) where this matrix enters the kinetic energy term.
The nonlinear susceptibility relates the secondorder induced density to the local electric field. The latter is a microscopic quantity that can be connected to the external field by knowing the linear response function. In the linear case it is a standard route to get the RPA dielectric response. In the nonlinear case the procedure is, probably, less known. Therefore, we follow here a very pedagogical treatment of Liebsch and Schaich Liebsch and Schaich (1989). In fact, they applied a trick suggested by Zangwill and Soven Zangwill and Soven (1980) to avoid the summation over the infinite number of unoccupied states for the calculation of the response functions. This allowed them to describe SHG at simple metal surfaces as effectively one dimensional systems (it is basically the same approach that enabled Prodan and Nordlander Prodan and Nordlander (2003) to treat very large spherical systems).
The outline of this work is as follows: In Sec. II we start with the most general relation between the induced densities and the effective fields [Eq. (2.5) of Liebsch and Schaich, Phys. Rev. B 40, 5401 (1989)] and formulate integral equations for the induced density with the noninteracting response functions as kernels. In Sec. III we discuss the case of a spherical symmetry and the simplifications it implies for the numerics. Finally, the second harmonic response of the magic Na cluster is presented in Sec. IV for the illustration of our theory. Based on our recently developed method for the solution of integral equations of this type we are able to drastically reduce the computational cost from to , where is the number of mesh points to represent the density.
Ii Integral equations
Our theory can easily be extended to include electron correlation effects by using the exchange correlation functional of DFT. Our formulation here is at the random phase approximation level. Within this approximation the densitydensity response function can be obtained as
(2) 
where is the Coulomb potential and is the noninteracting densitydensity response function (known as Lindhard function for the homogeneous electron gas in three dimensions, 3D):
(3) 
where is the Fermi function and , refer to the collections of quantum numbers that uniquely characterize the electronic states of the system. The infinitesimally small positive number shifts the poles from the real axis and ensures, thus, the causality of the response function. In what follows we will assume it can be incorporated in the variable.
Let us consider the response of the system subject to the harmonic electric field oscillating at the frequency , i.e. . Then, the induced density which oscillates at the frequency of the applied field is given by:
(4) 
where is the induced local field oscillating at the fundamental frequency and consisting of the external potential plus the Hartree potential corresponding to the induced density:
(5) 
The induced density oscillating at the double frequency results from the nonlinear process described by the response function and from the linear response to the local field oscillating at :
(6) 
Because there is no external field at the local field at this frequency is given by the Hartree potential:
(7) 
eq:n1,eq:v1 yield the integral equation for the linear density:
(8) 
while \crefeq:n2,eq:v2 result in the integral equation for the second harmonic density:
(9) 
We defined the source terms following Liebsch and Schaich (1989) as:
(10)  
(11) 
Eqs. (8, 10) and eqs. (9, 11) when coupled with an appropriate approximation for the noninteracting response functions allow to completely describe the linear and the secondharmonic response.
The semiclassical approximation for the has already been derived previously Pavlyukh et al. (2012) with the result:
(12) 
In Appendix A we derive the secondharmonic generation source term (49) that can be represented as:
(13) 
Although these expressions look rather complicated they can be further simplified in the case of spherical symmetry.
Iii Spherical symmetry
For spherically symmetric systems (i. e. ) \crefeq:i1,eq:i2 simplify considerably when projected onto spherical harmonics. Since these equations have the same functional form we introduce an index to label corresponding quantities. We use the multipole expansion of the fields:
This is a general form consistent with the spherical symmetry. Similar expansions will be used for the densities and the source terms:
For numerical calculations in this work we restrict ourselves to the mixture of the dipole and quadrupole excitations ( for ). Further simplifications can be achieved by considering the external field to be a planewave:
(14) 
and when we align the coordinate system along the direction:
(15) 
When the wavenumber is small compared to the dimension of the system (i. e. ) it is justified to replace the spherical bessel functions with their small argument approximations , use the definition of the multipole moments:
and express the external potential as:
(16)  
(17) 
where each term is a harmonic function i. e. . Thus, in this approximation, it is sufficient to consider only the first term in (12):
(18) 
With the help of \crefeq:n1,eq:v1 we have:
where we introduced the local plasmon frequency in analogy with the expression for homogenous systems . Finally, the secondorder source term can be given as a sum of two contributions :
(19a)  
(19b) 
The physical meaning of these two terms can be inferred by introducing the induced electric field in e.g. Eq. (13). The first term contains a contribution from the electric quadrupole term and from the density variation (the surface dipole term) , whereas the second term upon the use of the vector identity and one of the Maxwell equations can be interpreted as the magnetic “bulk” term and a surface term. According to the analysis of Wang, Cheng and Bower Wang et al. (1973) and Apell Apell (1983) the contribution from is small in the case of the nonlinear response from the surfaces and the incident electric field polarized in the plane of incidence. These arguments loose their power in the case of the nonlinear response from spherical objects. Thus, both terms must be considered. We will show below that the treatment of the second term is much more involved. Therefore, to illustrate our approach only the first term will be included into the numerical algorithm.
iii.1 Computational scheme
With these results in hands the numerical algorithm can be formulated as follows:

This integral equation is solved similarly to (8) in order to obtain the second order density.
iii.2 Linear response
In this subsection we focus on the first point of the program. Already on this level our approach is valuable as it provides the optical absorption crosssection.
The linear source term:
The integral kernel:
The response functions are invariant under the rotations of the system as a whole. Thus, in the linear case it can be expanded as:
Thus, the integral kernel of \crefeq:i1,eq:i2 can be projected onto the angular momentum eigenstates:
(21) 
It can be integrated using the spherical harmonics expression of the Coulomb potential [see Sec. 3.6 of Jackson (Ref. Jackson, 1999)]:
where () is the smaller (larger) of and . The gradient of the spherical harmonics need not to be considered in view of the symmetry of the density, the angular integration can be done beforehand and we obtain for (21):
(22) 
with the spherical pole Green function defined as
(23) 
With the help of (22) the integral equations (8) and (9) can be written in the unified form:
(24) 
This equation is central for our theory in both the linear and the nonlinear cases. An efficient method of its solution exists Pavlyukh et al. (2012). In the linear case the equation takes a more symmetric form when reformulated in terms of the polarizability function.
Observables:
The pole frequencydependent polarizability defined as:
can be computed as the response to the field . Let us introduce the positiondependent polarizability as:
and use the explicit form (20) for the linear source. Substituting these definitions in (24) we obtain the following integral equation:
(25) 
where
In the case of the dipolar response our result (24) coincides with Eq. (5) of Prodan and Nordlander (Ref. Prodan and Nordlander, 2003). Finally the frequencydependent polarizability is represented as the integral:
(26) 
iii.3 The nonlinear source term
The source term for the secondorder response is considerably more complicated. The central quantity is the induced potential . It can be evaluated by the integration of (5):
(27) 
where the density is obtained from the solution of (24). It is easy to see that in accordance with the earlier assumption it is sufficient to restrict ourselves to the cases of and . The most difficult part of the derivation: the transition from to can probably be obtained in closed form, however, for practical applications it is sufficient to have shorter form small solutions. They can be obtained with the maple computer algebra package:
(28) 
(29) 
where we introduced for brevity , , , .
The second part is presented in Appendix B for reference, however, will not be included into the numerical algorithm.
Iv Numerical results
In the present contribution we continue to study the optical properties of the magic Na cluster (Fig. 1a). This system simultaneously possesses completely closed geometric and electronic shells. This makes it exceptionally stable and attractive for the numerical calculations: its electronic structure can be obtained easily by using the density functional approach. We do not pursue here a fully atomistic approach, but rather solve the radial KohnSham equation (using the renormalized Numerov method) in the presence of the spherically symmetric ionic density. We either use the standard jellium model which starts from the homogeneous ionic density (Fig. 1b) or the density is obtained from the realistic geometric model by applying the gaussian broadening of Å width to ionic positions and performing a spherical averaging (Fig. 1c). In accordance with the proposed numerical scheme we present here a) the dipolar and the quadrupole linear optical absorption coefficients and the induced densities () (Fig. 2); b) results for the secondharmonic source ; and c) the nonlinear dipolar and quadrupole optical responses at frequency (Fig. 3).
While the linear optical response of metallic clusters is well understood: typically the optical absorption profile only slightly deviates from that of the idealistic system: a sphere filled with an electron gas of constant density. There, the optical absorption coefficient is peaked at the energies of the surface plasmon modes:
The position of the maxima is only weakly dependent on the details of the electronic structure: our calculations almost perfectly match corresponding idealistic values eV and eV for the bulk Na density (). The spilloff of the electron density in realistic systems mostly leads to the broadening of the surface plasmon resonances. To illustrate this fact we choose the offresonance value of the frequency ( eV) and plot the source and the induced density (Fig. 2). While in the idealistic situation the induced density is distributed over the surface of the sphere (where the electron density abruptly changes) in the realistic case we observe numerous features associated with slow electronic density variations within the cluster. They contribute to the optical absorption in the offresonance regime. However, the relative weight of these oscillations decreases when the frequency approaches the resonance. There, the fast density variation at the surface dominates the spectrum.
The nonlinear optical properties (Fig. 3) of even such simple systems are not fully understood. It is interesting to observe that two different excitation mechanisms (the quadrupole transition at or frequency) lead to almost identical frequency dependence. Unlike in the linear case, the efficiency of the frequency conversion vanishes at the plasmon resonance and has two pronounced peaks at the energy slightly below and above. We also do not observe a strong correlation between the source and the induced density as in the offresonant linear case (cf. blue and green curves). However, the spatial variation of these quantities is not erratic (as the plots might be suggesting). The complicated radial dependence is the result of the derivatives of the induced density and the potential in the nonlinear source terms. Thanks to the linearscaling algorithms we are using at each stage of the computation we are able to eliminate any spurious contributions from the numerical differentiation while ensuring the convergence of our results with respect to the number of discretization points and the value of the broadening parameter .
V Conclusions
We developed a semiclassical theory of the secondharmonic generation in spherical particles. Although it originates from the exact quantummechanical sumoverstates expression and takes the local fields into account it is free from the small system size limitation. Because an efficient method for the solution of the corresponding integral equations exists the theory can be applied even to macroscopic systems provided the electronic density can be found. For this purpose the jellium approximation for the ionic density can be used as we illustrate by computing the SHG spectrum of the Na cluster. The scattering crosssection and the angular distribution of the intensity are other important experimental observables. Calculation of these quantities will be presented elsewhere together with the inclusion of higher multipole moments allowing, thus, for the treatment of larger systems. The work along these lines is already in progress.
We believe that our method is sufficiently versatile as to allow for further extensions. In particular, it is straightforward to modify the method for systems with axial symmetry, or even to consider the symmetryfree case. To treat larger systems one must include higher multipole moments. This poses a question of how to find the secondorder source term in this case. We believe that a direct manipulation with the maple computer algebra system rather than a formal solution in terms of symbols is feasible.
Regarding the physical aspects of our approach: The fact that it is free of any adjustable parameters is not necessarily beneficial. In fact, for metallic systems with localized electrons peculiarities of the electronic structure might be reflected in the optical properties. In this case the classical electrodynamics approach with the experimentally measured dielectric function and susceptibility tensor might give more accurate results. On the other hand, for systems with simpler electronic structure our method is capable of taking into account quantum effects such as the spillout of the electronic density.
Finally, our work establishes the equivalence between the semiclassical approximation and the hydrodynamic approach of Apell Apell (1983). The latter, however, is just a classical theory if not corrected for local field effects. We also touched upon the question of the representability of the response functions in terms of electronic densities and show that in general a more complicated quantity such as the oneparticle density matrix must be introduced. However, for the secondharmonic generation these terms cancel in the final expression.
Acknowledgements.
The work is supported by DFGSFB762 and DFGSFB/TR 88. W. H. gratefully acknowledges the hospitality of the ”Nonequilibrium ManyBody Systems” group at MartinLutherUniversity HalleWittenberg during his sabbatical.Appendix A Semiclassical expression for the first hyperpolarizability
We start with the expression for the second order nonlinear response:
(30) 
where is the Fermi function and , , refer to collections of quantum numbers that uniquely characterize electronic states of the system. In what follows we use the the following notations , , etc. The infinitesimally small positive number shifts the poles from the real axis and assures, thus, the causality of the response function. In what follows we will assume that it can be incorporated in the variable.
Let us consider a generic type of integrals:
(31) 
and introduce our basic approximation. The semiclassical approximation (SCA) implies a high frequency condition . Thus, the above expression can be expanded as a power series of . The term proportional to in the expansion of vanishes in view of the completeness of the electron wavefunctions:
(32) 
In what follows we will assume all wavefunctions to be real. This can always be done without the loss of generality for timeinvariant systems. Odd power terms vanish in view of the symmetry of the expression with respect to permutation and .
The term has the following form:
(33) 
The first term vanishes after using the property (32), the integration over , and exploiting the permutation symmetry and of the expression under the integral. The term proportional to can be rewritten as and combined with the second terms. Finally our expression can be written as:
(34) 
with following notations:
In what follows we will make use of
(35) 
where . This follows from the Schrödinger equation.
Let us consider first the term. After the summation over and integration over we arrive at:
We further pull out of the integration and use the Gauss theorem to apply to the :
(36) 
Now the summation over and can be performed by introducing the oneparticle density matrix:
(37) 
with the diagonal elements given by the density:
(38) 
Quite general we can represent the density matrix in the form:
From the definition (37) a useful integration formula can be derived:
(39) 
Let us now return to (36) and perform the sum over the states using \crefeq:1n,eq:n.
(40) 
In order to evaluate the integral we must remove all differential operators from the function.
(41) 
It can be further simplified by using Green’s first identity for the first and third terms:
(42) 
Let us introduce the tensor and integrate:
(43) 
Using (35) we rewrite in the form without explicit energies:
(44) 
We perform the summations over the states and use the Gauss theorem to transfer on to obtain
(45) 
The integration yields: