A Semi-classical expression for the first hyperpolarizability

# Semi-classical approximation for the second harmonic generation in nanoparticles

## Abstract

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 Na cluster.

###### pacs:
78.67.-n, 73.20.Mf, 71.10.-w, 42.65.Ky, 31.15.A-,73.22.Dj

## 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 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 non-linear regimes Roke (2009). For electrons as a perturbation we refer to the overview  García de Abajo (2010). 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  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 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 Östling et al. (1993) and Dewitz, Hübner, and Bennemann Dewitz et al. (1996) proposed the anharmonic oscillator model to describe the second-order non-linear response from spherical systems. A further step was undertaken by Dadap et al. Dadap et al. (1999) who considered the small-particle 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 non-linear 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 non-linear Mie theory valid for particles of arbitrary sizes, de Beer and Roke included the sum-frequency 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 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 Prodan and Nordlander (2003). 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 semi-classical approach becomes very accurate. This is a marked finding as it allows to replace the complicated sum-over-states 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 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 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 sum-over-states expressions Pavlyukh et al. (2012) and using the assumption that the frequency of the external field is large compared with the single-particle 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:

 j(r,2ω)=\upalpha[f(E⋅∇)E+\upbeta(∇f)E2+\upgammaf(∇⋅E)E], (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 Rayleigh-Gans scattering approximation for a lattice of such particles Martorell et al. (1997), second harmonic generation by two-dimensional 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 sum-over-state quantum mechanical formula for the non-linear optical susceptibility Orr and Ward (1971) 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 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 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 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 one-particle density matrix whose diagonal elements are given by . This comes not as a surprise if we consider the analogy with the orbital-free kinetic density functional theory Wang and Carter (2002) 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 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 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 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.

We use atomic units, i. e., 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

 χ(r,r′;ω)=χ(0)(r,r′;ω)+∫dr1∫dr2χ(0)(r,r1;ω)×v(r1−r2)χ(r2,r′;ω), (2)

where is the Coulomb potential and is the non-interacting density-density response function (known as Lindhard function for the homogeneous electron gas in three dimensions, 3D):

 χ(0)(r,r′;ω)=2∑i,jfi−fjω+Ei−Ej+iη×ψi(r)ψ∗j(r)ψj(r′)ψ∗i(r′), (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:

 Extra open brace or missing close brace (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:

 φ(1)(r)=φ(0)(r)+∫dr′v(r−r′)δn(1)(r′). (5)

The induced density oscillating at the double frequency results from the non-linear process described by the response function and from the linear response to the local field oscillating at :

 δn(2)(r)=∫dr′∫dr′′χ(0)2(r;r′,r′′;ω)φ(1)(r′)φ(1)(r′′)+∫dr′χ(0)(r,r′;ω)φ(2)(r′). (6)

Because there is no external field at the local field at this frequency is given by the Hartree potential:

 φ(2)(r)=∫dr′v(r−r′)δn(2)(r′). (7)
\Cref

eq:n1,eq:v1 yield the integral equation for the linear density:

 δn(1)(r)=ξ(1)(r)+∫dr′∫dr′′χ(0)(r,r′;ω)v(r′−r′′)δn(1)(r′′), (8)

while \crefeq:n2,eq:v2 result in the integral equation for the second harmonic density:

 δn(2)(r)=ξ(2)(r)+∫dr′∫dr′′χ(0)(r,r′;ω)v(r′−r′′)δn(2)(r′′). (9)

We defined the source terms following Liebsch and Schaich (1989) as:

 ξ(1)(r) = ∫dr′χ(0)(r,r′;ω)φ(0)(r′), (10) ξ(2)(r) = ∫dr′∫dr′′χ(0)2(r;r′,r′′;ω)φ(1)(r′)φ(1)(r′′). (11)

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 second-harmonic response.

The semi-classical approximation for the has already been derived previously Pavlyukh et al. (2012) with the result:

 ξ(1)(r)=1ω2(∇n(r)⋅∇φ(0)(r)−n(r)Δφ(0)(r)). (12)

In Appendix A we derive the second-harmonic generation source term (49) that can be represented as:

 ξ(2)(r)=12ω4∇⋅(∇φ(1)(r)[n(r)Δφ(1)(r)+(∇φ(1)(r)⋅∇n(r))]+14n(r)∇(∇φ(1)(r))2). (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:

 φ(i)(r)=∑ℓmφ(i)ℓm(r)Yℓm(rr).

This is a general form consistent with the spherical symmetry. Similar expansions will be used for the densities and the source terms:

 δn(i)(r)=∑ℓmδn(i)ℓm(r)Yℓm(rr),ξ(i)(r)=1ω2i∑ℓmξ(i)ℓm(r)Yℓm(rr).

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 plane-wave:

 exp(ikr)=4π∑ℓmiℓjℓ(kr)Yℓm(Ωr)Yℓm(Ωk), (14)

and when we align the coordinate system along the -direction:

 exp(ikr)≡eikrcosθ=4π∑ℓiℓ√2ℓ+14πjℓ(kr)Yℓ0(rr). (15)

When the wave-number 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:

 Qℓm(r)=√4π2ℓ+1rℓYℓm(rr),

and express the external potential as:

 φ(0)(r) = ∑ℓ(ik)ℓ(2ℓ−1)!!Qℓ0(r),or (16) φ(0)ℓm(r) = (ik)ℓ(2ℓ−1)!!√4π2ℓ+1rℓδm0, (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):

 ξ(1)(r)=1ω2∇n(r)⋅∇φ(0)(r). (18)

With the help of \crefeq:n1,eq:v1 we have:

 nΔφ(1)=−ω2pδn(1),(∇φ(1)⋅∇n)=(ω2−ω2p)δn(1),

where we introduced the local plasmon frequency in analogy with the expression for homogenous systems . Finally, the second-order source term can be given as a sum of two contributions :

 ξ(2a)(r) = 12ω4∇⋅(δn(1)(r)(ω2−2ω2p(r))∇φ(1)(r)), (19a) ξ(2b)(r) = 18ω4∇⋅(n(r)∇(∇φ(1)(r))2). (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 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.

### iii.1 Computational scheme

With these results in hands the numerical algorithm can be formulated as follows:

• The integral equation (8) is solved with a source term  (18). It yields the first order density and, therefore, the local potential ;

• Eq. (19) is applied to generate the source term for the equation (9);

• 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 cross-section.

#### The linear source term:

Using (18) and the explicit form of the potential (17) we obtain:

 ξ(1)ℓm(r)=ℓn′(r)φ(0)ℓm(r)r, (20)

where denotes the derivative of the equilibrium density function with respect to the radial coordinate.

#### 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:

 χ(r,r′;ω)=∑lμχlμ(r,r′;ω)Y∗lμ(Ω)Ylμ(Ω′).

Thus, the integral kernel of \crefeq:i1,eq:i2 can be projected onto the angular momentum eigenstates:

 Kℓ(r,r′)=∫dΩY∗ℓm(Ω)×∫dΩ′∇n(r)⋅∇rv(|r−r′|)Yℓm(Ω′). (21)

It can be integrated using the spherical harmonics expression of the Coulomb potential [see Sec. 3.6 of Jackson (Ref. Jackson, 1999)]:

 v(r−r′)=1|r−r′|=∑ll∑μ=−l4π2l+1rlY∗lμ(Ω)Ylμ(Ω′),

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):

 Kℓ(r,r′)=4π2ℓ+1n′(r)∂∂r(rℓ)=−4π2ℓ+1n′(r)rℓ−1r′ℓ+1Gsphℓ(r,r′), (22)

with the spherical -pole Green function defined as

 Gsphℓ(r,r′)=(ℓ+1)θ(r−r′)(r′r)2ℓ+1−ℓθ(r′−r). (23)

With the help of (22) the integral equations (8) and (9) can be written in the unified form:

 (1−ω2p(r)ω2)δn(i)ℓm(r;ω)=ξ(i)ℓ(r;ω)ω2i−4π2ℓ+1n′(r)ω2rℓ−1×∫∞0dr′(1r′)ℓ−1Gsphℓ(r,r′)δn(i)ℓm(r′;ω). (24)

This equation is central for our theory in both the linear and the non-linear cases. An efficient method of its solution exists Pavlyukh et al. (2012). In the linear case the equation takes a more symmetric form when re-formulated in terms of the polarizability function.

#### Observables:

The -pole frequency-dependent polarizability defined as:

 αℓm(ω)=−∫dr∫dr′Q∗ℓm(r)χ(r,r′;ω)Qℓm(r′),

can be computed as the response to the field . Let us introduce the position-dependent polarizability as:

 αℓ(r;ω)=√4π2ℓ+1δn(1)ℓ0(r;ω)

and use the explicit form (20) for the linear source. Substituting these definitions in (24) we obtain the following integral equation:

 αℓ(r;ω)=α(0)ℓ(r;ω)×[ℓ−∫∞0dr′(1r′)ℓ−1Gsphℓ(r,r′)αℓ(r′;ω)]. (25)

where

 α(0)ℓ(r;ω)=−4π2ℓ+1rℓ−1n′(r)ω2−ω2p(r).

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 frequency-dependent polarizability is represented as the integral:

 αℓ0(ω)=∫∞0drrℓ+2α(0)ℓ(r;ω). (26)

### iii.3 The non-linear source term

The source term for the second-order response is considerably more complicated. The central quantity is the induced potential . It can be evaluated by the integration of (5):

 φ(1)ℓm(r)=φ(0)ℓm(r)+4π2ℓ+1∫dr′r′2rℓδn(1)ℓm(r′), (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:

 ξ(2a)1=12√5π[u′(n1f′2+n2f′1)+3ur2(n1f2+n2f1)+u(n′1f′2+n′2f′1−8πn1n2)], (28)
 Missing or unrecognized delimiter for \Big (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 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  Å 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 second-harmonic source ; and c) the non-linear 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:

 ωℓ=√ℓ2ℓ+1ωp.

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 spill-off 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 ( 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 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 spectrum.

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 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 discretization points and the value of the broadening parameter .

## V Conclusions

We developed a semi-classical theory of the second-harmonic 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 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 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 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 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 spill-out of the electronic density.

Finally, our work establishes the equivalence between the semi-classical 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 one-particle density matrix must be introduced. However, for the second-harmonic generation these terms cancel in the final expression.

###### Acknowledgements.
The work is supported by DFG-SFB762 and DFG-SFB/TR 88. W. H. gratefully acknowledges the hospitality of the ”Nonequilibrium Many-Body Systems” group at Martin-Luther-University Halle-Wittenberg during his sabbatical.

## Appendix A Semi-classical expression for the first hyperpolarizability

We start with the expression for the second order non-linear response:

 χ(0)2(r;r′,r′′;ω)=∑i,j,kψ∗k(r)ψi(r)ψ∗i(r′)ψj(r′)ψ∗j(r′′)ψk(r′′)2ω+Ei−Ek+iη×(fk−fjω+Ej−Ek+iη−fj−fiω+Ei−Ej+iη) (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:

 Φ(r)=∫d(r′r′′)χ(0)2(r;r′,r′′;ω)ϕ(r′)ϕ(r′′) (31)

and introduce our basic approximation. The semi-classical 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 wave-functions:

 ∑iψ∗i(r1)ψi(r2)=δ(r1−r2). (32)

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 and .

The term has the following form:

 Φ(4)(r)=∫d(r′r′′)∑i,j,kXi,j,k(r,r′,r′′)ϕ(r′)ϕ(r′′)×[fkjE2kj−12EikfkjEkj+14E2ikfkj]. (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 re-written as and combined with the second terms. Finally our expression can be written as:

 Φ(4)(r)=14Φ(4a)(r)−Φ(4b)(r)+12Φ(4c)(r) (34)

with following notations:

 Φ(4a)(r) = ∫d(r′r′′)∑i,j,kXi,j,k(r,r′,r′′)[E2ikfk]ϕ(r′)ϕ(r′′), Φ(4b)(r) = ∫d(r′r′′)∑i,j,kXi,j,k(r,r′,r′′)[EikEjkfj]ϕ(r′)ϕ(r′′), Φ(4c)(r) = ∫d(r′r′′)∑i,j,kXi,j,k(r,r′,r′′)[EikEjkfk]ϕ(r′)ϕ(r′′).

In what follows we will make use of

 Eijψi(r)ψj(r)=−12∇⋅ξij(r), (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:

 Φ(4a)(r)=14∫dr′ϕ2(r′)∑i,kfk[∇r⋅ξik(r)][∇r′⋅ξik(r′)].

We further pull out of the integration and use the Gauss theorem to apply to the :

 Φ(4a)=−∇r4⋅∫dr′∑i,kfkξik(r)[ξik(r′)⋅∇r′ϕ2(r′)]. (36)

Now the summation over and can be performed by introducing the one-particle density matrix:

 ∑ifiψi(r1)ψi(r2)=γ(r1,r2), (37)

with the diagonal elements given by the density:

 ∑ifiψi(r)ψi(r)=n(r), (38)

Quite general we can represent the density matrix in the form:

 γ(r1,r2)=f(r1)f(r2)g(|r1−r2|).

From the definition (37) a useful integration formula can be derived:

 ∫dr′δ(r−r′)F(r′)∇rγ(r,r′)=∫dr′δ(r−r′)F(r′)∇r′γ(r,r′)=12F(r)∇rn(r). (39)

Let us now return to (36) and perform the sum over the states using \crefeq:1n,eq:n.

 Φ(4a)=−14∑ab∇ar∫dr′∇r′ϕ2(r′)×[γ(r,r′)∇ar∇br′δ(r,r′)−∇arγ(r,r′)∇br′δ(r,r′)−∇arδ(r,r′)∇br′γ(r,r′)+δ(r,r′)∇ar∇br′γ(r,r′)]. (40)

In order to evaluate the integral we must remove all differential operators from the -function.

 Extra open brace or missing close brace (41)

It can be further simplified by using Green’s first identity for the first and third terms:

 Φ(4a)=Δ4∫dr′δ(r,r′)[Δϕ2(r′)γ(r,r′)+2∇br′ϕ2(r′)∇br′γ(r,r′)]−12∑a∇ar∫dr′δ(r,r′)[Δ[ϕ2(r′)]∇arγ(r,r′)+2∑b∇ar∇br′γ(r,r′)∇br′ϕ2(r′)]. (42)

Let us introduce the tensor and integrate:

 Φ(4a)=14Δ[nΔϕ2+(∇ϕ2⋅∇n)]−14∇⋅[∇nΔϕ2]−∑ab∇a[Tab∇bϕ2]. (43)

Using (35) we re-write in the form without explicit energies:

 Φ(4b)(r)=∇r4⋅∫d(r′r′′)ϕ(r′)ϕ(r′′)∑i,j,kfjξik(r)[∇r′′⋅ξjk]×ψi(r′)ψj(r′). (44)

We perform the summations over the states and use the Gauss theorem to transfer on to obtain

 Extra open brace or missing close brace (45)

The integration yields:

 Φ(4b1)(r) = −14∑ab∇ar[∇brϕ(r)∇ar∫dr′ϕ(r′)δ(r,r′)∇brγ(r,r′) +14∑ab∇ar[∇brϕ(r)∫dr′ϕ(r′)δ(r,r′)∇ar∇brγ(r,r′)] = −14∑a∇a[12(∇n⋅∇ϕ)∇aϕ+ϕ∑bTab∇bϕ], Φ(4b3)(r) = 18∑a∇a[ϕ∇a(∇ϕ⋅∇n)−2ϕ∑bTab∇bϕ], Φ(4b4)(r) = −14∑ab∇ar[ϕ(r)∫dr′′∇br′′ϕ(r′′)γ(r′′,r)∇br′′∇arδ(r,r′′)] = −14∑ab∇ar[ϕ(r)∇ar∫dr′′∇br′′ϕ(r′′)γ(r′′,r)∇br′′δ(r,r′′)] + 14∑ab∇ar[ϕ(r)∫dr′′∇br′′ϕ