Semi-classical approximation for the second harmonic generation in nanoparticles

Second harmonic generation by spherical nanoparticles is a non-local optical process that can also be viewed as the result of the non-linear response of the thin interface layer. The classical electrodynamic description, based e.g. on the non-linear Mie theory, entails the knowledge of the dielectric function and the surface non-linear optical susceptibility, both quantities are usually assumed to be predetermined, for instance from experiment. We propose here an approach based on the semi-classical approximation for the quantum sum-over-states expression that allows to capture the second-order optical process from first principles. A key input is the electronic density, which can be obtained from effective single particle approaches such as the density-functional 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 Na2869 cluster.


I. 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 k to the system [1]. To name but few examples, we mention here the response of nanoparticles to light in the linear [2] and the non-linear regimes [3]. For electrons as a perturbation we refer to the overview [4]. A well-studied route to address these issues is the concept of the dielectric response functions which is detailed in standard books such as in [1], but also in recent reviews, e.g. in Refs. [2][3][4].
For extended bulk metals and metallic surfaces the dielectric response was extensively studied (Ref. [5,6] and references therein). An extension to finite-size 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 [7] and Dewitz, Hübner, and Bennemann [8] proposed the anharmonic oscillator model to describe the second-order nonlinear response from spherical systems. A further step was undertaken by Dadap et al. [9] who considered the small-particle limit (ka < 1) and described SHG as a mixture of dipole and quadrupole excitation processes (a 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 χ (2) (ω). In this approach the inversion symmetry is not a necessary assumption as it is * yaroslav.pavlyukh@physik.uni-halle.de 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 [10] developed the non-linear Mie theory valid for particles of arbitrary sizes, de Beer and Roke included the sum-frequency generation mechanism into the considerations [11], the cylindrical geometry was treated by Dadap [12] and finally the theory for arbitrarily shaped particles was developed by de Beer, Roke, and Dadap [13].
All these theories based on classical electrodynamics rely on the knowledge of the frequency-dependent dielectric function and the non-linear 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 semi-infinite and finite size systems on the nanometer range. To shed light on these issues, it is desirable to have a quantum theory for the non-linear 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 [14]. They succeeded to push the limits of the time-dependent 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 sum-over-states quantum mechanical expression [15] 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 [16,17] for fullerenes. The usefulness of the jellium model was demonstrated by the pioneering works of Ekardt on sodium clusters [18,19] or of Puska and Nieminen on C 60 [20].) we can obtain the ground state electronic density from the solution of the Kohn-Sham equations and express the response function in its terms. The Drude dielectric function follows automatically.
The semi-classical approximation is rooted in the work of Mukhopadhyay and Lundqvist [21] 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 60 [22] or carbon nanotubes [23]. The equation can also be derived starting from the quantum mechanical sum-overstates expressions [24] and using the assumption that the frequency of the external field is large compared with the singleparticle gap (ω E g ) [25].
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 [26] classically treated second harmonic generation from alkali metals. A decade later Apell [27] derived an expression for the second harmonic current in the form: where for the unperturbed ground state electron density in the form n(r) = n 0 f (r) the α parameter is a function of n 0 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 [28], the Rayleigh-Gans scattering approximation for a lattice of such particles [29], second harmonic generation by two-dimensional particles [30], or a more recent treatment of arbitrary geometries [31]. 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 (E) 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 sum-over-state quantum mechanical formula for the non-linear optical susceptibility [15] and employing the high-frequency 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 high-frequency semi-classical expansion. It is interesting to notice that the 1/ω 4 asymptotic behavior of this quantity is at variance with the results of Scan-dolo and Bassani [32] who predict a 1/ω 6 decay. This seems, however, to be a direct consequence of the assumption of the all-dipole excitation mechanism that underly their study. It is possible to prove from general principles that the inclusion of the quadrupole excitations leads to the 1/ω 4 asymptotics [33].
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 [34] SHG comprises three processes in which the 2ω 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 n(r). As we demonstrate below, the answer is negative, one additionally needs the one-particle density matrix γ(r, r ) whose diagonal elements are given by n(r). This comes not as a surprise if we consider the analogy with the orbital-free kinetic density functional theory [35] where this matrix enters the kinetic energy term.
The non-linear susceptibility relates the second-order 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 non-linear case the procedure is, probably, less known. Therefore, we follow here a very pedagogical treatment of Liebsch and Schaich [36]. In fact, they applied a trick suggested by Zangwill and Soven [37] 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 [14] 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 non-interacting 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 − 2869 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 O(N 3 ) to O(N), where N is the number of mesh points to represent the density.
We use atomic units, i. e., = e = m e = 4πε 0 = 1 throughout. Two appendices contain mathematical details to make the exposition in Secs. II-IV self-contained.

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 density-density response function can be obtained as where v(r−r 1 ) is the Coulomb potential and χ (0) (r, r ; ω) is the non-interacting density-density response function (known as Lindhard function for the homogeneous electron gas in three dimensions, 3D): where f is the Fermi function and i, j 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. ϕ (0) (r; t) = ϕ (0) (r) e −iωt . Then, the induced density which oscillates at the frequency of the applied field is given by: (4) where ϕ (1) (r) 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: The induced density oscillating at the double frequency results from the non-linear process described by the χ (0) 2 response function and from the linear response to the local field ϕ (2) (r) oscillating at 2ω: Because there is no external field at 2ω the local field at this frequency is given by the Hartree potential: Equations (4) and (5) yield the integral equation for the linear density: while eqs. (6) and (7) result in the integral equation for the second harmonic density: We defined the source terms following [36] as: Eqs. (8,10) and eqs. (9,11) when coupled with an appropriate approximation for the non-interacting response functions allow to completely describe the linear and the secondharmonic response. The semi-classical approximation for the ξ (1) (r) has already been derived previously [24] with the result: In Appendix A we derive the second-harmonic generation source term (A20) that can be represented as: 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. n(r) ≡ n(r)) eqs. (8) and (9) simplify considerably when projected onto spherical harmonics. Since these equations have the same functional form we introduce an index i = 0, 1, 2 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 (ϕ (0) m = 0 for > 2). Further simplifications can be achieved by considering the external field to be a plane-wave: and when we align the coordinate system along the kdirection: When the wave-number k is small compared to the dimension of the system a (i. e. ka < 1) it is justified to replace the spherical bessel functions with their small argument approximations j (z) ∼ z n /(2 + 1)!!, use the definition of the multipole moments: and express the external potential as: where each term is a harmonic function i. e. ∆Q m (r) = 0. Thus, in this approximation, it is sufficient to consider only the first term in (12): With the help of eqs. (4) and (5) we have: where we introduced the local plasmon frequency in analogy with the expression for homogenous systems ω 2 p (r) = 4πn(r). Finally, the second-order source term can be given as a sum of two contributions ξ (2) (r) = ξ (2a) (r) + ξ (2b) (r): The physical meaning of these two terms can be inferred by introducing the induced electric field E = −∇ϕ (1) in e.g. Eq. (13). The first term contains a contribution from the electric quadrupole term E(∇ · E) and from the density variation (the surface dipole term) E(E · ∇n), whereas the second term upon the use of the vector identity 1 2 ∇E 2 = E×∇×E+(E·∇)E and one of the Maxwell equations E × ∇ × E = iω/cE × H can be interpreted as the magnetic "bulk" term and a surface term. According to the analysis of Wang, Cheng and Bower [26] and Apell [27] the contribution from ξ (2b) is small in the case of the non-linear 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 non-linear 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.

A. Computational scheme
With these results in hands the numerical algorithm can be formulated as follows: i) The integral equation (8) is solved with a source term (18). It yields the first order density δn (1) (r) and, therefore, the local potential ϕ (1) (r); ii) Eq. (19) is applied to generate the source term for the equation (9); iii) This integral equation is solved similarly to (8) in order to obtain the second order density.

B. 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 cross-section.
a. The linear source term: Using (18) and the explicit form of the potential (17) we obtain: where n (r) denotes the derivative of the equilibrium density function with respect to the radial coordinate. b. 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 eqs. (8) and (9) can be projected onto the angular momentum eigenstates: where r < (r > ) is the smaller (larger) of r and r . 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): with the spherical -pole Green function defined as With the help of (22) the integral equations (8) and (9) can be written in the unified form: This equation is central for our theory in both the linear and the non-linear cases. An efficient method of its solution exists [24]. In the linear case the equation takes a more symmetric form when re-formulated in terms of the polarizability function. c. Observables: The -pole frequency-dependent polarizability defined as: can be computed as the response to the field φ (0) (r) = 4π 2 +1 r . Let us introduce the position-dependent polarizability as: and use the explicit form (20) for the linear source. Substituting these definitions in (24) we obtain the following integral equation: where .
In the case of the dipolar response our result (24) coincides with Eq. (5) of Prodan and Nordlander (Ref. 14). Finally the frequency-dependent polarizability is represented as the integral: C. The non-linear source term The source term for the second-order response is considerably more complicated. The central quantity is the induced potential ϕ (1) (r). It can be evaluated by the integration of (5): 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 = 1, 2 and m = 0. The most difficult part of the derivation: the transition from ϕ (1) (r) to ξ (2) (r) 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: where we introduced for brevity . The second part ξ (2b) (r) 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 − 2869 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 Kohn-Sham 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 ∼ 1 Å 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 α 0 (ω) and the induced densities δn (1) 0 (r; ω) ( = 1, 2) (Fig. 2); b) results for the second-harmonic source ξ (2) (r); and c) the non-linear dipolar and quadrupole optical responses at 2ω 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 ω 1 = 3.45 eV and ω 2 = 3.78 eV for the bulk Na density (r s = 3.96). 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 off-resonance value of the frequency (ω * = 5.75 eV) and plot the source ξ (2) (r, ω * ) and the induced density δn (2) (r, ω * ) (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 off-resonance 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 spec-trum.
The non-linear 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 2ω 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 off-resonant 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 non-linear source terms. Thanks to the linear-scaling 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 dis-  cluster for the two models at Fig. 2. Right panels show the spatial dependence of the corresponding source terms ξ (2) (r, ω * ) and the resulting induced densities δn (2) (r, ω * ) for a particular frequency value ω * = 5.75 eV. cretization points and the value of the broadening parameter η.

V. CONCLUSIONS
We developed a semi-classical theory of the secondharmonic generation in spherical particles. Although it originates from the exact quantum-mechanical sum-over-states 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 − 2869 cluster. The scattering cross-section 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 mo-ments 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 symmetry-free case. To treat larger systems one must include higher multipole moments. This poses a question of how to find the second-order 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 3 j 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 d-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 spill-out of the electronic density.
Finally, our work establishes the equivalence between the semi-classical approximation and the hydrodynamic approach of Apell [27]. 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 one-particle density matrix must be introduced. However, for the second-harmonic generation these terms cancel in the final expression. 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: Φ(r) = d(r r )χ (0) 2 (r; r , r ; ω) φ(r ) φ(r ) (A2) and introduce our basic approximation. The semi-classical approximation (SCA) implies a high frequency condition |E i − E j | ω. Thus, the above expression can be expanded as a power series of ω. The term proportional to 1/ω 2 in the expansion of Φ(r) vanishes in view of the completeness of the electron wave-functions: i ψ * i (r 1 )ψ i (r 2 ) = δ(r 1 − r 2 ).
In what follows we will assume all wave-functions to be real. This can always be done without the loss of generality for time-invariant systems. Odd power terms vanish in view of the symmetry of the expression with respect to permutation i k and r r . The 1/ω 4 term has the following form: The first term vanishes after using the property (A3), the integration over r , and exploiting the permutation symmetry j k and r r of the expression under the integral. The term proportional to E 2 ik f j can be re-written as 2E ik E jk f j and combined with the second terms. Finally our expression can be written as: with following notations: Φ (4a) (r) = d(r r ) i, j,k X i, j,k (r, r , r ) E 2 ik f k φ(r ) φ(r ), i, j,k X i, j,k (r, r , r ) E ik E jk f j φ(r ) φ(r ), In what follows we will make use of