text:
# Formalism for Multiphoton Plasmon Excitation in Jellium Clusters.

## 1 Introduction

In the present paper we demonstrate that, in addition to dipole plasmon excitations multipole plasmons (quadrupole, octupole, etc) contribute to the multiphoton excitation process, which results in a significant difference of plasmon resonance profiles between the cross sections for multiphoton and single-photon absorption. We have developed a formalism from which the cross sections for multiphoton excitation can be worked out. The balance between multipole surface and volume plasmon contributions is analyzed. Our results are obtained within a theoretical model for the multiphoton excitation of a jellium cluster. This model is applicable to metal clusters, to fullerenes and to any type of cluster in which a strong delocalization of valence atomic orbitals takes place. The theoretical formalism we have developed is not confined in its application to photons. It can also be used to describe any kind of higher order plasmon excitation processes, for example those which arise by multiple scattering of electrons within the cluster.

Recently, a number of papers have discussed metallic clusters  and fullerenes  in strong laser fields. The theoretical approach usually followed is to solve a time-dependent local density equation numerically (LDA or local-density approximation), and the regime most commonly studied involves intense, short laser pulses, for which the turn-on and turn-off properties significantly affect excitation.

We report here on a different problem: our initial interest lies in lower laser powers, for which multiphoton excitation just begins to intrude, and we are interested in developing a formalism to describe the interaction between collective modes and a laser field in the multiphoton regime. For this purpose, a semiclassical model, in which the collective flow of charge is driven by a periodic field, is established, and we relate it to the multiphoton absorption cross section of the cluster, which takes account of quantum mechanics. Of course, one can in principle extend this treatment to consider the turn-on and turn-off of laser pulses, to treat the interaction numerically for various power levels and initial charge distributions. What we wish to point out, however, is a novel feature, which arises even for an infinite wavetrain interacting with a cluster (the simplest and most fundamental problem): multiple plasmon excitations are driven by multiphoton excitations. In the present paper, we explain by what mechanism this arises.

Our approach is based on different principles from the LDA. Instead of using the Kohn-Sham formalism, we propose a hydrodynamic approach. In the LDA, all quantities are made to depend solely on charge density, and currents are subsequently made to appear by solving a time-dependent equation. In our approach, we start out from the continuity equation and the Euler equation, from which both the current flow and the density are obtained in a completely consistent way. Our theory is local, and does not include exchange. In principle, it would be possible to extend it by building in exchange and correlations in a similar way to the LDA. One of the benefits of our approach is that all the momentum transfer terms are included in the formulation, which leads to the presence of both volume and surface plasmon terms. If one wishes to simplify the theory, it is possible towards the end of the calculation to assume zero momentum transfer, in which case the volume plasmon terms disappear from the problem.

Surface plasmon excitations are well known in atomic cluster physics. The dipole surface plasmons are responsible for the formation of giant resonances in photoabsorption spectra of metal clusters (see e.g. ). They also play an important role in inelastic collisions of charged particles with metal clusters . The role of surface plasmon excitations in inelastic electron-cluster scattering was thoroughly studied in , and it was demonstrated that collective excitations make a significant contribution to the electron energy loss spectrum (EELS) in the region of the surface plasmon resonance. With increasing scattering angle, plasmon excitations of higher angular momenta become more and more prominent.

Plasmons are characteristic of delocalised electrons, and therefore the jellium model provides the most appropriate starting point for a discussion of plasmon excitation in the multiphoton regime. This regime is particularly appropriate for the study of atomic clusters, which are rather fragile objects, and readily explode under very strong irradiation

The inelastic scattering of fast electrons on metal clusters in the range of transferred energies above the ionization threshold was considered in . It was demonstrated that, in this energy range, volume plasmons dominate the contribution to the differential cross section, resulting in a resonance behaviour. The volume plasmon resonances excited in the cluster during a collision decay via the ionization process. The resonance frequency and the autoionization width of the volume plasmon excitations have both been determined in .

The role of the polarization interaction and plasmon excitations in the process of electron attachment to metal clusters has also been examined both theoretically  and experimentally . It was demonstrated that plasmon excitations induce a resonance enhancement of the electron attachment cross section.

Our paper is organized as follows. In section 2, we derive quantum mechanical expressions for the cross sections of multipole (quadrupole, octupole etC) plasmon excitations taking place in the multiphoton absorption regime and estimate the cross sections on the basis of the plasmon resonance approximation. We present and discuss the multiphoton absorption profiles and demonstrate a significant change of the multiphoton absorption profiles in the cross sections for multiphoton as compared to single-photon absorption. In section 3, we establish certain connections of the cross sections with the variation of electron density in the cluster due to the external field and present the formalism for the calculation of the electron density variation in the cluster due to the external field, based on the use of the hydrodynamic Euler equation and on the equation of continuity. We apply the general formalism to the description of fast electron-cluster scattering and multiphoton absorption. In section 4, we calculate the multipole moments of the system induced by the external field on the basis of the formalism outlined in section 3. We analyse the plasmon resonance structure of the induced multipole moments and conclude that it is analogous to the one arising in the multiphoton absorption cross sections calculated in section 2. In section 5, we draw conclusions from this work. In A, we show how matrix elements for collective transitions can be calculated on the basis of the sum rules. In B, we present details of calculations of the angular integrals arising in the formalism outlined here.

## 2 Plasmon resonance approximation for multiphoton absorption cross sections

First, we consider the cross section for multiphoton absorption in jellium clusters and demonstrate that multipole plasmon excitations are essential to this process. The discussion in this section is based on the plasmon resonance approximation, which is introduced below.

### 2.1 Single photon absorption

Let us start by considering the simplest example and calculate the cross section for single-photon absorption in the plasmon resonance approximation.

The single-photon absorption cross section in the dipole approximation reads as

$$\sigma _1=\frac{4\pi ^2e^2}{c}\omega \underset{n}{}|z_{on}|^2\delta (\omega _{no}\mathrm{}\omega )$$
(1)
Here, $`e`$ is the charge of electron, $`c`$ is the velocity of light, $`\mathrm{}`$ is Planck’s constant, $`\omega _{no}=\epsilon _n\epsilon _0`$ is the electron excitation energy, $`\omega `$ is the photon frequency and $`ez_{on}`$ is the matrix element of the $`z`$\- component of the cluster dipole moment. The summation over $`n`$ includes all final states of the excited electron, which belong to both the discrete and the continuous spectra.

In the jellium picture, which works reasonably well for metal clusters and to some extent for fullerenes, the main contribution to the cross section (1) arises from a small group of excited states or sometimes even from a single transition of frequency close to the classical Mie resonance frequency – also known as the frequency of the plasmon resonance. For a spherical metal cluster, this frequency is given by (see e.g.  and section 4 of this paper)

$$\omega _l^2=\frac{4\pi Ne^2}{mV}\frac{l}{(2l+1)}.$$
(2)
Here $`V=4\pi R^3/3`$ is the cluster volume, where $`R=r_oN^{1/3}`$ is the cluster radius, $`r_o`$ is the Wigner-Seitz radius; $`N`$ is the number of delocalized electrons in a cluster, $`l`$ is the angular momentum of the plasmon mode, $`m`$ is the electron mass. Note that, by using a single photon of energy $`14eV`$, one can, in practice, excite only $`l=1`$ dipole plasmon oscillations in a metallic cluster.

For nearly-spherical fullerenes $`C_{20}`$ or $`C_{60}`$, the plasmon resonance frequency is equal to (see )

$$\omega _l^2=\frac{l(l+1)N}{(2l+1)R^3}$$
(3)
where $`N`$ is the total number of delocalised electrons (4 electrons per atom times the number of carbon atoms in the fullerene molecule) and $`R`$ is the the radius of the fullerene.

The plasmon resonance approximation is based on the fact that excitations in the vicinity of a plasmon resonance exhaust the sum rule almost completely (see ), which means that the summation in the sum rule (see e.g. )

$$\underset{n}{}\omega _{no}|z_{on}|^2=\frac{N\mathrm{}^2}{2m}$$
(4)
need be performed only over excited states in the vicinity of the Mie resonance.

Now, assuming a Lorentzian distribution of width $`\mathrm{\Gamma }_1`$ for the plasmon resonance states and replacing the delta function, $`\delta (\omega _{no}\mathrm{}\omega )`$, in (1) by the profile (see e.g. )

$$\delta (\omega _{no}\mathrm{}\omega )\frac{\mathrm{\Gamma }_1}{2\pi \mathrm{}((\omega _1\omega )^2+\frac{\mathrm{\Gamma }_1^2}{4})},$$
(5)
one recovers the well-known expression for the single-photon absorption cross section (see e.g. )

$$\sigma _1=\frac{\pi Ne^2}{mc}\frac{\mathrm{\Gamma }_1}{(\omega _1\omega )^2+\frac{\mathrm{\Gamma }_1^2}{4}}\frac{4\pi Ne^2}{mc}\frac{\omega ^2\mathrm{\Gamma }_1}{(\omega _1^2\omega ^2)^2+\omega ^2\mathrm{\Gamma }_1^2}$$
(6)
The width $`\mathrm{\Gamma }_1`$ is due to Landau damping. Its calculation for metal clusters is performed, for example, in .

The cross section (6) reproduces correctly the appearence of the plasmon resonances in single-photon absorption spectra of metal clusters and fullerenes, although some details of the experimentally observed profiles are naturally beyond the plasmon resonance approximation. The discussion of these details is not the scope of the present work. This can only be done accurately enough on the basis of ab initio many-body theories. Instead, we analyse the multiphoton absorption cross sections on the basis of the plasmon resonance approximation and elucidate the role of multipole plasmon excitations in its formation, because our interest lie in establishing the physical mechanisms which underpin multiphoton excitation.

### 2.2 Two-photon absorption

In the dipole approximation, the two-photon absorption cross section is equal to

$$\sigma _2=\frac{32\pi ^3e^4\mathrm{}}{c^2}\omega ^2\underset{n}{}\left|\underset{m}{}\frac{z_{nm}z_{mo}}{\mathrm{}\omega \omega _{mo}+\mathrm{i}\delta }\right|^2\delta (\omega _{no}2\mathrm{}\omega )$$
(7)
We evaluate the cross section (7) in the same way as for the single-photon case. The main contribution to the sum over the intermediate states $`m`$ arises from the virtual dipole plasmon excitations. Therefore, one derives

$$\underset{n}{}\left|\underset{m}{}\frac{z_{nm}z_{mo}}{\mathrm{}\omega \omega _{mo}+\mathrm{i}\delta }\right|^2\frac{N}{2m\mathrm{}\omega _1}\underset{n}{}\frac{|z_{n1}|^2}{(\omega \omega _1)^2+\mathrm{\Gamma }_1^2/4}$$
(8)
Here, we have also introduced a dipole plasmon resonance width $`\mathrm{\Gamma }_1`$ and used the sum rule (4) for the evaluation of the matrix elements for the dipole plasmon excitation $`|r_{10}|^2N\mathrm{}/2m\omega _1`$. The remaining matrix elements $`z_{n1}`$ in (8) describe dipole transitions from the dipole plasmon resonance state to other excited states. Matrix elements for these transitions obey the dipole selection rule. This means that the angular momentum of the final state can only be equal to either $`l=0`$ or $`l=2`$. According to (2) and (3), there is no surface plasmon excitation with $`l=0`$ either in metal clusters or in fullerenes. Thus, only transitions to the states with $`l=2`$ are of interest.

These arguments show that, by using two photons simultaneously, one can excite the quadrupole plasmon resonance in a metal cluster or in a fullerene with a frequency given in (2) and (3) respectively. When calculating the cross section (7) in the vicinity of the quadrupole plasmon resonance excitation, i.e. at $`2\omega \omega _2`$, it is sufficient to consider only transitions to the resonance final state, i.e. to put $`_n|z_{n1}|^2|z_{21}|^2`$ (here and below we use indices 1 and 2 to designate the dipole and quadrupole plasmon resonance states) and to replace the delta function, $`\delta (\omega _{no}2\mathrm{}\omega )`$ by a Lorentzian distribution of width $`\mathrm{\Gamma }_2`$ (see e.g. )

$$\delta (\omega _{no}2\mathrm{}\omega )\frac{\mathrm{\Gamma }_2}{2\pi \mathrm{}((\omega _22\omega )^2+\frac{\mathrm{\Gamma }_2^2}{4})}.$$
(9)
By substituting (8) and (9) in (7), one derives

$`\sigma _2=\left({\displaystyle \frac{4\pi Ne^2}{mc}}\right)^2{\displaystyle \frac{m|z_{21}|^2}{2N\mathrm{}\omega _1}}{\displaystyle \frac{\omega ^2}{(\omega \omega _1)^2+\frac{\mathrm{\Gamma }_1^2}{4}}}{\displaystyle \frac{\mathrm{\Gamma }_2}{(\omega _22\omega )^2+\frac{\mathrm{\Gamma }_2^2}{4}}}`$ (10)
This result demonstrates that the photoabsorption profile in the two-photon process differs substantially from the single-photon case. The cross section (10) has a resonance at the dipole plasmon frequency and, in addition, also contains the quadrupole plasmon resonance at $`\omega =\omega _2/2`$.

The cross section (10) depends on $`|z_{21}|^2`$. The transition matrix element $`z_{21}`$ describes the electron transition between the dipole and quadrupole plasmon resonance states. This is a single electron transition rather than a collective one. Therefore, calculation of $`z_{21}`$ on the basis of the sum rule (4) would lead to a significant overestimate of the value of this matrix element. Instead, one can use Heisenberg’s uncertainty principle for the evaluation of $`z_{21}`$ . By estimating the radial component of the momentum of a single electron in a dipole and quadrupole plasmon oscillatory mode as $`p_1m\omega _1\mathrm{\Delta }R`$ and $`p_2m\omega _2\mathrm{\Delta }R`$ respectively, one derives

$$|z_{21}|min(\mathrm{\Delta }z_1,\mathrm{\Delta }z_2)\frac{\mathrm{}}{p_2}A\frac{\mathrm{}}{m\omega _1\mathrm{\Delta }R}$$
(11)
Here, $`A`$ is a dimensionless constant, of the order of one, $`\mathrm{\Delta }z_1\mathrm{}/p_1`$ and $`\mathrm{\Delta }z_2\mathrm{}/p_2`$ are the uncertainties relating to an electron in the dipole and quadrupole plasmon modes respectively, $`\mathrm{\Delta }R`$ is the width of the domain in the vicinity of the cluster surface within which plasmon excitations take place. In A we prove the correctness of this estimate and demonstrate that matrix element $`z_{21}`$ is equal to

$$z_{21}=\frac{8}{3}\left(\frac{6}{5}\right)^{1/4}\frac{\mathrm{}}{m\omega _1\mathrm{\Delta }R}$$
(12)
By substituting (12) in (10), we obtain the final expression for the two-photon absorption cross section:

$`\sigma _2=\left({\displaystyle \frac{4\pi Ne^2}{mc}}\right)^2{\displaystyle \frac{A^2\mathrm{}}{2m\omega _1N\mathrm{\Delta }R^2}}{\displaystyle \frac{\omega ^2}{\omega _1^2}}{\displaystyle \frac{1}{(\omega \omega _1)^2+\frac{\mathrm{\Gamma }_1^2}{4}}}{\displaystyle \frac{\mathrm{\Gamma }_2}{(\omega _22\omega )^2+\frac{\mathrm{\Gamma }_2^2}{4}}}`$ (13)
where $`A=\frac{8}{3}\left(\frac{6}{5}\right)^{1/4}2.79`$.

We note that the cross section (13) depends explicitly on Planck’s constant $`\mathrm{}`$, while the cross section (6) does not. The independence of (6) from $`\mathrm{}`$ is connected with the fact that plasmon oscillations are a purely classical effect, while the dependence of (13) on $`\mathrm{}`$ arises from the interaction between dipole and quadrupole plasmon modes as can be seen from estimate (11) and the explicit expression (12). This indicates that it is meaningful to treat plasmon excitations classically, while the coupling of various plasmon modes in the multiphoton photoabsorption process must be treated beyond purely classical theory.

In figure 1 we plot the cross section profiles per unit atom for single-photon (dashed line) and two-photon (solid line) absorption calculated according to (6) and (13). These profiles do not depend on the number of atoms in the cluster. Note that the scales are not the same for the two curves for reasons of definition of the cross sections in the single- and two-photon cases, but are both given in atomic units. The peak in the single-photon plot gives the location of the dipole resonance. The other peak in the two-photon plot is the quadrupole resonance. This figure demonstrates a significant difference between the nature of the profiles, arising from the presence of quadrupole plasmon excitation in the two-photon case. In this calculation we have input $`r_0=4.0`$ and $`\mathrm{\Gamma }_1=\omega _1/4`$, $`\mathrm{\Gamma }_2=\omega _2/4`$, $`\mathrm{\Delta }R=r_0`$. The choice of these parameters can be different for different clusters, but it should always lead to qualitatively similar single- and two-photon absorption profiles. An accurate determination of the parameters is only possible on the basis of ab initio calculations.

### 2.3 $`n`$-photon absorption

The formalism we have developed can also be used for the evaluation of the multiphoton absorption cross sections for a larger number of photons. In the dipole approximation the n-photon absorption cross section has the following form:

$$\sigma _n=\frac{(2\pi )^{n+1}n!^2e^{2n}\mathrm{}^{n1}}{c^n}\omega ^n\underset{k}{}\left|M_k\right|^2\delta (\omega _{ko}n\mathrm{}\omega )$$
(14)
Here the amplitude $`M_n`$ is equal to

$`M_k`$ $`=`$ $`{\displaystyle \underset{m_{n1}}{}}{\displaystyle \underset{m_{n2}}{}}\mathrm{}.{\displaystyle \underset{m_1}{}}{\displaystyle \frac{z_{km_{n1}}}{((n1)\mathrm{}\omega \omega _{m_{n1}o}+\mathrm{i}\delta )}}`$ (15)
$``$ $`{\displaystyle \frac{z_{m_{n1}m_{n2}}}{((n2)\mathrm{}\omega \omega _{m_{n2}o}+\mathrm{i}\delta )}}.\mathrm{}{\displaystyle \frac{z_{m_10}}{(\mathrm{}\omega \omega _{m_1o}+\mathrm{i}\delta )}}`$
The plasmon resonance structure of the multiphoton absorption cross section (15) can be analysed in a similar way to the previous treatments for the single- and two-photon cases. This analysis immediately leads to the important conclusion that plasmon resonances with larger angular momenta (octupole, etc) can be excited in the multiphoton absorption regime. Thus, for example, with three photons, the octupole plasmon resonance at $`\omega =\omega _3/3`$ will also be excited. This analysis, however, leaves undefined the matrix elements for electronic transitions between various plasmon modes. Estimates of these can be performed either on the basis of Heisenberg’s Uncertainty Principle or by a calculation similar to the one for $`Z_{12}`$ (see (12) and (10)) in the two-photon case, but their accurate evaluation is not trivial.

Note that the plasmon resonance approximation allows one to analyse only the plasmon resonance excitations that are characterised by relatively low angular momenta, because electron excitations in the cluster with large angular momenta $`l`$ have single-particle character. This follows, for instance, from the fact that with increasing $`l`$ the wave length of the surface plasmon mode, $`2\pi R/l`$, becomes smaller than the characteristic wave length of the delocalised electrons at the Fermi surface, $`2\pi \mathrm{}/\sqrt{2m\mathrm{\Delta }\epsilon }`$, where $`\mathrm{\Delta }\epsilon `$ is the characteristic electron excitation energy in the cluster. In other words, excitations with angular momenta comparable with the characteristic electron angular momenta of the ground state exhibit single-particle rather than collective character. Therefore, when analyzing contributions of the plasmon resonance modes to the multiphoton absorption cross section, one should consider only the lowest angular momenta. For example according to the jellium model, the maximum angular momentum of the delocalised electrons in the $`Na_{40}`$ cluster is equal to 4. Therefore, only dipole and quadrupole collective modes can be expected in this case. With increasing cluster size the number of essential plasmon modes grows as $`R`$.

## 3 Hydrodynamic description of collective motion of the electron density in a cluster

The multipole plasmon resonances arising in the multiphoton absorption cross sections, should also appear in other physical characteristics of the cluster, which can be probed in the multiphoton absorption regime. In the situation where plasmon resonance excitations are the dominant contribution to the multiphoton absorption cross section, it is natural to seek and analyse the plasmon resonance structure of the variation of electron density induced by the radiation field. The variation of electron density is a characteristic of the system, allowing one easily to connect classical and quantum descriptions of the excitation process, because charge density variation has the same meaning in both quantum and classical mechanics. A classical description of the electron density variation in a cluster is appropriate in the situation where plasmon excitations dominate over the single-particle spectrum, because plasmon oscillations in clusters are an essentially classical effect.

### 3.1 Basic equations

Since our object of interest is the excitation of plasmons in metallic clusters, which have a distinctly classical nature, we now describe the collective motion of the electron density using Euler’s equation and the equation of continuity.

Euler’s equation couples the acceleration of the electron density $`d𝐯/dt`$ with the total local electric field $`𝐄`$ acting on the density at the point $`(𝐫,t)`$. It has the following form:

$$\frac{d𝐯(𝐫,t)}{dt}=\frac{e}{m}𝐄(𝐫,t)$$
(16)
The electric field $`𝐄`$ includes both the external field acting on the cluster and the polarization contribution arising from the variation of electron density. Expressing the total derivative on the left hand side of (16) as the sum of two contributions, arising from the change in velocity of the electron density in time and in space, one obtains:

$$\frac{𝐯(𝐫,t)}{t}+\{𝐯(𝐫,t)\}𝐯(𝐫,t)=\frac{e}{m}\phi (𝐫,t)\frac{e}{m}𝑑𝐫^{}\frac{\delta \rho (𝐫^{},t)}{|𝐫𝐫^{}|}$$
(17)
Here $`\phi (𝐫,t)`$ is the potential of the external field. The second term on the right hand side of (17) describes the polarization force due to the variation of electron density $`\delta \rho (𝐫,t)`$.

We assume that the external potential $`\phi (𝐫,t)`$ is the solution of the wave equation. Therefore, we can put

$$\phi (𝐫,t)=\mathrm{e}^{i\omega t}\phi (𝐫),$$
(18)
where $`\phi (𝐫,t)`$ satisfies the equation

$$\mathrm{\Delta }\phi (𝐫)=k^2\phi (𝐫).$$
(19)
Here $`k=\omega /c`$, $`c`$ is the velocity of light, but in principle one can postulate a more complex dispersion law. We need consider only the positive frequency solution of the wave equation, because the formalism for the negative frequency solution is analogous to it.

The total electron density in the cluster is equal to

$$\rho (𝐫,t)=\rho _o(𝐫)+\delta \rho (𝐫,t),$$
(20)
where $`\rho _o(𝐫)`$ is the electron density distribution in a free cluster without an external field and $`\delta \rho (𝐫,t)`$ is the variation of electron density caused by the external field and the polarization force acting together.

The motion of electron density in the cluster obeys the equation of continuity, which reads:

$$\frac{\rho (𝐫,t)}{t}+\{\rho (𝐫,t)𝐯(𝐫,t)\}=0$$
(21)
The simultaneous solution of equations (17), (20) and (21) with appropriate initial conditions and the initial distribution $`\rho _o(𝐫)`$ allow one to determine the variation of electron density $`\delta \rho (𝐫,t)`$ as well as its velocity $`𝐯(𝐫,t)`$. We solve this problem by using a perturbative approach on the external field $`\phi (𝐫,t)`$.

### 3.2 Perturbation theory

It is easy to estimate the relative value of the first and the second terms on the left hand side of (17). We see that the second term is negligible, provided the condition $`Em\omega ^2R/e`$ is fulfilled. Substituting here the characteristic values $`\omega 0.1`$, $`R10`$, one derives $`E0.1`$ in atomic units or $`E510^8B/cm`$. Below, we assume that this condition is fulfilled and neglect the second term on the left hand side of (17), which means physically that the external field causes only a small spatial inhomogeneity in the electron density distribution within the cluster. In this limit, Euler’s equation reduces to a Newtonian equation, which describes electronic motion in the cluster under the action of the external field and the polarization force.

We express the solutions of (17) and (21) in the following form:

$$\delta \rho (𝐫,t)=\underset{n=1}{\overset{\mathrm{}}{}}\delta \rho _n(𝐫)\mathrm{e}^{in\omega t}$$
(22)
$$𝐯(𝐫,t)=\underset{n=1}{\overset{\mathrm{}}{}}𝐯_n(𝐫)\mathrm{e}^{in\omega t}$$
(23)
By substituting these expansions into (17), (21) and performing simple transformations, one derives

$$𝐯_n(𝐫)=\frac{ie}{mn\omega }\delta _{n,1}\phi (𝐫)+\frac{ie}{mn\omega }𝑑𝐫^{}\frac{\delta \rho _n(𝐫^{})}{|𝐫𝐫^{}|}$$
(24)
$$i\omega n\delta \rho _n(𝐫)+\left\{\rho _0(𝐫)𝐯_n(𝐫)\right\}+\underset{k^{}=1}{\overset{n1}{}}\left\{\delta \rho _k^{}(𝐫)𝐯_{nk^{}}(𝐫)\right\}=0$$
(25)
Here $`\delta _{n,1}`$ is the Kroneker symbol. One can exclude $`𝐯_n(𝐫)`$ from equation (25) by the substitution of (24) in (25). Performing this transformation with the simultaneous use of (19) and $`\mathrm{\Delta }|𝐫𝐫^{}|^1=4\pi \delta (𝐫𝐫^{})`$, one derives the following equation:

$`\left((\omega n)^2{\displaystyle \frac{4\pi e}{m}}\rho _o(𝐫)\right)\delta \rho _n(𝐫)+{\displaystyle \frac{e}{m}}\rho _o(𝐫){\displaystyle 𝑑𝐫^{}\frac{\delta \rho _n(𝐫^{})}{|𝐫𝐫^{}|}}=`$
$`={\displaystyle \frac{e}{m}}\delta _{n,1}\left(\rho _o(𝐫)\phi (𝐫)k^2\phi (𝐫)\rho _o(𝐫)\right)+`$
$`+i\omega {\displaystyle \underset{k^{}=1}{\overset{n1}{}}}(\delta \rho _k^{}(𝐫)𝐯_{nk^{}}(𝐫))`$ (26)
The left hand side of equation (26) describes eigen-oscillations of the electron density. The electron density is almost constant within the cluster but varies rapidly near the cluster surface. Therefore, the terms proportional to $`\rho _o(𝐫)`$ and $`\rho _o(𝐫)`$ on the left hand side of (26) determine the square of the frequency of the volume and surface plasmon oscillations respectively. The right hand side in (26) describes a driving force acting on the eigen- plasmon oscillations.

The set of non-linear equations (24) and (26) must be solved iteratively. It is clear from the form of the equations that the index $`n`$ corresponds to the order of perturbation theory on the external field $`\phi (𝐫)`$.

Indeed, for $`n=1`$ equation (26) reduces to

$`\left(\omega ^2{\displaystyle \frac{4\pi e}{m}}\rho _o(𝐫)\right)\delta \rho _1(𝐫)+{\displaystyle \frac{e}{m}}\rho _o(𝐫){\displaystyle 𝑑𝐫^{}\frac{\delta \rho _1(𝐫^{})}{|𝐫𝐫^{}|}}=`$
$`={\displaystyle \frac{e}{m}}\left(\rho _o(𝐫)\phi (𝐫)k^2\phi (𝐫)\rho _o(𝐫)\right)`$ (27)
and, for $`n=2`$, one derives

$`\left((2\omega )^2{\displaystyle \frac{4\pi e}{m}}\rho _o(𝐫)\right)\delta \rho _2(𝐫)+{\displaystyle \frac{e}{m}}\rho _o(𝐫){\displaystyle 𝑑𝐫^{}\frac{\delta \rho _2(𝐫^{})}{|𝐫𝐫^{}|}}=`$
$`=i\omega (\delta \rho _1(𝐫)𝐯_1(𝐫))`$ (28)
Equations (27) and (28) show that the variation $`\rho _1(𝐫^{})`$ describes the linear response of the electron subsystem to the the external field $`\phi (𝐫)`$, while $`\rho _2(𝐫^{})`$ arises only in the second order of perturbation theory on $`\phi (𝐫)`$, because $`\rho _1(𝐫^{})\phi (𝐫)`$ and $`𝐯_1(𝐫)\phi (𝐫)`$.

Solving the set of equations (24) and (26) with $`\phi (𝐫)`$ describing the dipole electron-photon interaction up to the $`n`$-th order, one can calculate the variation of electron density in the cluster caused by the field of $`n`$ photons.

The set of equations (24) and (26) is not confined in its application to photons. It can also be used to describe the dynamics of electron density under the action of any kind of external field, for example the electric field of a charged projectile colliding with the cluster. Indeed, by considering the partial spherical harmonic of the Fourier image of the Coulomb field of the projectile particle, one can derive from (27) the same expression for the variation of the electron density $`\delta \rho _1(𝐫)`$ as follows from the purely electrodynamical perturbative approach to the electron scattering problem .

### 3.3 Spherically symmetric case

Equations (24) and (26) are valid for an arbitrary shape of the initial distribution $`\rho _o(𝐫)`$. In the spherical case, the angular parts in (24) and (26) can be separated. Thus, the cross section for $`n`$-photon absorption can in principle be extracted for arbitrarily large n, although the calculations become more and more tedious the higher $`n`$ is. Let us consider this formalism in more detail.

In the case of the spherically symmetric initial distribution, one can put $`\rho _o(𝐫)=\rho _0(r)`$. This relationship allows one easily to exclude angular variables from equation (26). Using this relationship together with the partial expansion for $`\delta \rho _n(𝐫)`$ and $`\phi (𝐫)`$,

$$\delta \rho _n(𝐫)=\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=l}{\overset{l}{}}\delta \rho _{l,m}^n(r)Y_{l,m}(𝐧_r)$$
(29)
$$\phi (𝐫)=\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=l}{\overset{l}{}}\phi _{l,m}(r)Y_{l,m}(𝐧_r),$$
(30)
one derives

$`\left((\omega n)^2{\displaystyle \frac{4\pi e}{m}}\rho _o(r)\right)\delta \rho _{l,m}^n(r)+{\displaystyle \frac{4\pi e\rho _o^{}(r)}{m(2l+1)}}{\displaystyle 𝑑r^{}G_l(r,r^{})\delta \rho _{l,m}^n(r^{})}=`$
$`={\displaystyle \frac{e}{m}}\delta _{n,1}(\rho _o(r)\phi _{l,m}(r)k^2\phi _{l,m}^{}(r)\rho _o^{}(r)){\displaystyle \frac{e}{m}}{\displaystyle \underset{j=1}{\overset{n1}{}}}{\displaystyle \frac{1}{nj}}{\displaystyle }d\mathrm{\Omega }_{𝐧_r}Y_{l,m}^{}(𝐧_r)\times `$
$`\times \{\delta \rho _j(𝐫)(\delta _{nj,1}\phi (𝐫)+{\displaystyle }d𝐫^{}{\displaystyle \frac{\delta \rho _{nj}(𝐫^{})}{|𝐫𝐫^{}|}})`$
$`\delta \rho _j(𝐫)(\delta _{nj,1}\phi (𝐫)k^2+4\pi \delta \rho _{nj}(𝐫))\}`$ (31)
Here, $`Y_{l,m}(𝐧_r)`$ is the sperical harmonic corresponding to the angular momentum $`l`$ and the projection of the angular momentum $`m`$. When deriving (31), we have multiplied both sides of equation (26) by the sperical harmonic $`Y_{l,m}^{}(𝐧_r)`$ and then integrated over $`d\mathrm{\Omega }_{𝐧_r}`$. We also used the well known expansion (see e.g. )

$$\frac{1}{|𝐫𝐫^{}|}=\frac{1}{2l+1}\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=l}{\overset{l}{}}B_l(r,r^{})Y_{l,m}(𝐧_r)Y_{l,m}^{}(𝐧_r^{})$$
(32)
where function $`B_l(r,r^{})`$ is defined as follows

$$B_l(r,r^{})=\frac{r^l}{r^{l+1}}\mathrm{\Theta }(r^{}r)+\frac{r^l}{r^{l+1}}\mathrm{\Theta }(rr^{})$$
(33)
Here $`\mathrm{\Theta }(r^{}r)`$ is the step function. In equation (31) we have introduced the function $`G_l(r,r^{})=r^2B_l(r,r^{})/r`$, which is of the the form

$$G_l(r,r^{})=l\frac{r^{l1}}{r^{l1}}\mathrm{\Theta }(r^{}r)(l+1)\frac{r^{l+2}}{r^{l+2}}\mathrm{\Theta }(rr^{})$$
(34)
When deriving (31), we have also made obvious transformations of the sum over $`k^{}`$ on the right hand side of equation (26) using (19) and (24). Note that the sum over $`j`$ in (31) still contains the integrals over the angular variables. The integration over the angular variables in the sum is straightforward, but somewhat cumbersome. It is also clear that the non-gradient terms in the sum contain the integration of the product of three spherical harmonics, if one expands $`\delta \rho _{nj}(𝐫)`$ and $`\phi (𝐫)`$ according to (29) and (30), which we denote as

$$I_1(l,m|l_1,m_1,|l_2,m_2)=d\mathrm{\Omega }_{𝐧_r}Y_{l,m}^{}(𝐧_r)Y_{l_1,m_1}(𝐧_r)Y_{l_2,m_2}(𝐧_r)$$
(35)
The gradient terms in the sum contain the integration of a spherical harmonic multiplied by the scalar product of the two vector spherical harmonics:

$`I_2(l,m|l_1,m_1,|l_2,m_2)=\sqrt{l_1(l_1+1)l_2(l_2+1)}{\displaystyle }d\mathrm{\Omega }_{𝐧_r}\times `$
$`Y_{l,m}^{}(𝐧_r)𝐘_{l_1,m_1}^{(1)}(𝐧_r)𝐘_{l_2,m_2}^{(1)}(𝐧_r).`$ (36)
This type of integral arises, when one expresses the gradients of the potential and the density according to (see e.g. )

$`\delta \rho _{l_1,m_1}^j(r)Y_{l_1,m_1}(𝐧_r)=`$
$`=\delta \rho _{l_1,m_1}^j(r)𝐘_{l_1,m_1}^{(1)}(𝐧_r)+\sqrt{l_1(l_1+1)}{\displaystyle \frac{1}{r}}\delta \rho _{l_1,m_1}^j(r)𝐘_{l_1,m_1}^{(1)}(𝐧_r)`$ (37)
$`\mathrm{\Phi }_{l_2,m_2}^{nj}(r)Y_{l_2,m_2}(𝐧_r)=`$
$`=\mathrm{\Phi }_{l_2,m_2}^{(nj)}(r)𝐘_{l_2,m_2}^{(1)}(𝐧_r)+\sqrt{l_2(l_2+1)}{\displaystyle \frac{1}{r}}\mathrm{\Phi }_{l_2,m_2}^{(nj)}(r)𝐘_{l_2,m_2}^{(1)}(𝐧_r)`$ (38)
Here $`𝐘_{l,m}^{(1)}(𝐧_r)`$, $`𝐘_{l,m}^{(1)}(𝐧_r)`$ are, respectively, the longitudinal and the transverse vector spherical harmonics, the definition of which one can find in . We mention some properties of these vector hamonics: $`𝐘_{l,m}^{(1)}(𝐧_r)=𝐧_𝐫Y_{l,m}(𝐧_r)`$, $`𝐘_{l,m}^{(1)}(𝐧_r)=_𝛀Y_{l,m}(𝐧_r)/\sqrt{l(l+1)}`$ and $`𝐘_{l,m}^{(1)}(𝐧_r)𝐘_{l,m}^{(1)}(𝐧_r)=0`$.

The potential $`\mathrm{\Phi }_{l_2,m_2}^{nj}(r)`$ in (38) is as follows

$$\mathrm{\Phi }_{l_2,m_2}^{nj}(r)=\delta _{nj,1}\phi _{l_2,m_2}(r)+\frac{4\pi }{(2l_2+1)}𝑑r^{}r^2B_l(r,r^{})\delta \rho _{l_2,m_2}^{nj}(r^{})$$
(39)
Since the explicit expressions for $`I_1(l,m|l_1,m_1,|l_2,m_2)`$ and $`I_2(l,m|l_1,m_1,|l_2,m_2)`$ are somewhat lengthy, they are presented in B.

Using the formulae written above, one can easily rewrite (31) in the following form

$`\left((\omega n)^2{\displaystyle \frac{4\pi e}{m}}\rho _o(r)\right)\delta \rho _{l,m}^n(r)+{\displaystyle \frac{4\pi e\rho _o^{}(r)}{m(2l+1)}}{\displaystyle 𝑑r^{}G_l(r,r^{})\delta \rho _{l,m}^n(r^{})}=`$
$`={\displaystyle \frac{e}{m}}\delta _{n,1}\left(\rho _o(r)\phi _{l,m}(r)k^2\phi _{l,m}^{}(r)\rho _o^{}(r)\right)`$
$`{\displaystyle \frac{e}{m}}{\displaystyle \underset{j=1}{\overset{n1}{}}}{\displaystyle \frac{1}{nj}}{\displaystyle \underset{\stackrel{l_1,m_1}{l_2,m_2}}{}}I_1(l,m|l_1,m_1,|l_2,m_2)\times `$
$`\times \left\{\delta \rho ^j(r)\mathrm{\Phi }_{l_2,m_2}^{(nj)}(r)\delta \rho _{l_1,m_1}^j(r)(\delta _{nj,1}\phi _{l_2,m_2}(r)k^2+4\pi \delta \rho _{l_2,m_2}^{nj}(r))\right\}`$
$`{\displaystyle \frac{e}{mr^2}}{\displaystyle \underset{j=1}{\overset{n1}{}}}{\displaystyle \frac{1}{nj}}{\displaystyle \underset{\stackrel{l_1,m_1}{l_2,m_2}}{}}I_2(l,m|l_1,m_1,|l_2,m_2)\delta \rho _{l_1,m_1}^j(r)\mathrm{\Phi }_{l_2,m_2}^{nj}(r)`$ (40)
Here and below, we assume that the summation over $`l_1,m_1`$ and $`l_2,m_2`$ is performed with the same limits as in (29) and (30).

### 3.4 Surface and volume plasmons

We now analyse equation (40) and demonstrate that it describes both surface and volume plasmon oscillations. The surface and volume solutions of (40) can be separated from each other, if one assumes that the initial distribution of electron density has the form

$$\rho _o(r)=\frac{Ne}{V}\mathrm{\Theta }(R_vr)$$
(41)
Here $`N`$ is the total number of delocalized electrons in the cluster volume $`V=4\pi R_v^3/3`$.

In this case it is natural to look for the solution of (40) by expressing it in the following form

$$\delta \rho _{l,m}^{(n)}(r)=\delta \rho _{l,m}^{s(n)}\delta (R_sr)+\delta \rho _{l,m}^{v(n)}(r)\mathrm{\Theta }(R_vr)$$
(42)
In (41) and (42) we have introduced the two radii $`R_v`$ and $`R_s`$ and assumed that $`R_v<R_s=R`$, but $`R_vR_s=R`$, where $`R`$ is the cluster radius. Such a relationship is necessary for the elimination of the uncertainties, which arise in (42) and the subsequent formulae in the vicinity of the cluster radius.

Substituting (41) and (42) in (40), performing straighforward but lengthy calculations of the integrals, one derives the following equation

$`\left((\omega n)^2\omega _p^2\right)\delta \rho _{l,m}^{v(n)}(r)\mathrm{\Theta }(Rr)+`$
$`+\left((\omega n)^2\omega _l^2\right)\delta \rho _{l,m}^{s(n)}\delta (Rr)=`$
$`={\displaystyle \frac{4\pi eN(l+1)}{m(2l+1)VR^{l+2}}}\delta (Rr){\displaystyle _0^R}𝑑r^{}r^{l+2}\delta \rho _{l,m}^{v(n)}(r^{})+`$
$`+{\displaystyle \frac{Ne^2}{mV}}\delta _{n,1}\phi _{l,m}(r)k^2\mathrm{\Theta }(Rr)+{\displaystyle \frac{Ne^2}{mV}}\delta _{n,1}\phi _{l,m}^{}(R)\delta (Rr)+`$
$`+{\displaystyle \frac{e}{m}}\delta (Rr){\displaystyle \underset{j=1}{\overset{n1}{}}}{\displaystyle \frac{1}{nj}}{\displaystyle \underset{\stackrel{l_1,m_1}{l_2,m_2}}{}}I_1(l,m|l_1,m_1,|l_2,m_2)\times `$
$`\times \{\delta \rho _{l_1,m_1}^{s(j)}(\delta _{nj,1}\phi _{l_2,m_2}^{\prime \prime }(R)+{\displaystyle \frac{4\pi }{(2l_2+1)}}{\displaystyle _0^R}dr^{}D_{l_2}(R,r^{})\delta \rho _{l_2,m_2}^{v(nj)}(r^{}))+`$
$`+\delta \rho _{l_1,m_1}^{v(j)}(R)\left(\delta _{nj,1}\phi _{l_2,m_2}^{}(R)+{\displaystyle \frac{4\pi }{(2l_2+1)}}{\displaystyle _0^R}𝑑r^{}G_{l_2}(R,r^{})\delta \rho _{l_2,m_2}^{v(nj)}(r^{})\right)+`$
$`+\delta _{nj,1}k^2\delta \rho _{l_1,m_1}^{s(j)}\phi _{l_2,m_2}(R)\}`$
$`{\displaystyle \frac{e}{m}}\mathrm{\Theta }(Rr){\displaystyle \underset{j=1}{\overset{n1}{}}}{\displaystyle \frac{1}{nj}}{\displaystyle \underset{\stackrel{l_1,m_1}{l_2,m_2}}{}}I_1(l,m|l_1,m_1,|l_2,m_2)\times `$
$`\times \{\delta \rho _{l_1,m_1}^{v(j)}(r)(\delta _{nj,1}\phi _{l_2,m_2}^{}(r)+{\displaystyle \frac{4\pi }{(2l_2+1)}}{\displaystyle _0^R}dr^{}G_{l_2}(r,r^{})\delta \rho _{l_2,m_2}^{v(nj)}(r^{}))`$
$`\delta _{nj,1}k^2\delta \rho _{l_2,m_2}^{v(j)}(r)\phi _{l_2,m_2}(r)4\pi \delta \rho _{l_2,m_2}^{v(j)}(r)\delta \rho _{l_2,m_2}^{v(j)}(r)\}`$
$`{\displaystyle \frac{e}{mR^2}}\delta (Rr){\displaystyle \underset{j=1}{\overset{n1}{}}}{\displaystyle \frac{1}{nj}}{\displaystyle \underset{\stackrel{l_1,m_1}{l_2,m_2}}{}}I_2(l,m|l_1,m_1,|l_2,m_2)\times `$
$`\times \delta \rho _{l_1,m_1}^{s(j)}\left(\delta _{nj,1}\phi _{l_2,m_2}(R)+{\displaystyle \frac{4\pi }{(2l_2+1)}}{\displaystyle _0^R}𝑑r^{}r^2B_{l_2}(R,r^{})\delta \rho _{l_2,m_2}^{v(nj)}(r^{})\right)`$
$`{\displaystyle \frac{e}{mr^2}}\mathrm{\Theta }(Rr){\displaystyle \underset{j=1}{\overset{n1}{}}}{\displaystyle \frac{1}{nj}}{\displaystyle \underset{\stackrel{l_1,m_1}{l_2,m_2}}{}}I_2(l,m|l_1,m_1,|l_2,m_2)\times `$
$`\times \delta \rho _{l_1,m_1}^{v(j)}(r)\left(\delta _{nj,1}\phi _{l_2,m_2}(r)+{\displaystyle \frac{4\pi }{(2l_2+1)}}{\displaystyle _0^R}𝑑r^{}r^2B_{l_2}(r,r^{})\delta \rho _{l_2,m_2}^{v(nj)}(r^{})\right)`$ (43)
When deriving (43), we have used the fact that $`\mathrm{\Theta }^{}(Rr)=\delta (Rr)`$. Also we have introduced the function

$$D_l(r,r^{})=l_2(l_21)\frac{r^{l_22}}{r^{l_21}}\mathrm{\Theta }(r^{}r)+(l_2+1)(l_2+2)\frac{r^{l_2+2}}{r^{l+3}}\mathrm{\Theta }(rr^{})$$
(44)
The left hand side of equation (43) describes volume and surface eigen-oscillations of the electron density characterised by the angular momentum $`l`$. The surface plasmon resonance frequency $`\omega _l`$ is the same as in (2). The volume plasmon resonance frequency is equal to

$$\omega _p=\sqrt{\frac{4\pi e^2N}{mV}}$$
(45)
In equation (43), $`\omega _p`$ appears in expressions involving $`l`$, but is nevertheless is independent of $`l`$, as one sees in equation (45). The physical reason for this is that the volume plasmon oscillation is degenerate with $`l`$. The right hand side of (43) describes a driving force acting on the eigen-plasmon oscillations.

Surface and volume terms on the right hand side of equation (43) have not been regrouped, in order to stress their correspondencce with terms in (40). It is seen from (43) that equations for the volume and surface plasmon oscillations can be separated and will then read as follows:

$`\left((\omega n)^2\omega _p^2\right)\delta \rho _{l,m}^{v(n)}(r)={\displaystyle \frac{Ne^2}{mV}}\delta _{n,1}\phi _{l,m}(r)k^2`$
$`{\displaystyle \frac{e}{m}}{\displaystyle \underset{j=1}{\overset{n1}{}}}{\displaystyle \frac{1}{nj}}{\displaystyle \underset{\stackrel{l_1,m_1}{l_2,m_2}}{}}I_1(l,m|l_1,m_1,|l_2,m_2)\times `$
$`\times \{\delta \rho _{l_1,m_1}^{v(j)}(r)(\delta _{nj,1}\phi _{l_2,m_2}^{}(r)+{\displaystyle \frac{4\pi }{(2l_2+1)}}{\displaystyle _0^R}dr^{}G_{l_2}(r,r^{})\delta \rho _{l_2,m_2}^{v(nj)}(r^{}))`$
$`\delta _{nj,1}k^2\delta \rho _{l_2,m_2}^{v(j)}(r)\phi _{l_2,m_2}(r)4\pi \delta \rho _{l_2,m_2}^{v(j)}(r)\delta \rho _{l_2,m_2}^{v(j)}(r)\}`$
$`{\displaystyle \frac{e}{mr^2}}{\displaystyle \underset{j=1}{\overset{n1}{}}}{\displaystyle \frac{1}{nj}}{\displaystyle \underset{\stackrel{l_1,m_1}{l_2,m_2}}{}}I_2(l,m|l_1,m_1,|l_2,m_2)\times `$
$`\times \delta \rho _{l_1,m_1}^{v(j)}(r)\left(\delta _{nj,1}\phi _{l_2,m_2}(r)+{\displaystyle \frac{4\pi }{(2l_2+1)}}{\displaystyle _0^R}𝑑r^{}r^2B_{l_2}(r,r^{})\delta \rho _{l_2,m_2}^{v(nj)}(r^{})\right)`$ (46)
$`\left((\omega n)^2\omega _l^2\right)\delta \rho _{l,m}^{s(n)}={\displaystyle \frac{Ne^2}{mV}}\delta _{n,1}\phi _{l,m}^{}(R)`$
$`{\displaystyle \frac{4\pi eN(l+1)}{m(2l+1)VR^{l+2}}}{\displaystyle _0^R}𝑑r^{}r^{l+2}\delta \rho _{l,m}^{v(n)}(r^{})+`$
$`+{\displaystyle \frac{e}{m}}{\displaystyle \underset{j=1}{\overset{n1}{}}}{\displaystyle \frac{1}{nj}}{\displaystyle \underset{\stackrel{l_1,m_1}{l_2,m_2}}{}}I_1(l,m|l_1,m_1,|l_2,m_2)\times `$
$`\times \{\delta \rho _{l_1,m_1}^{s(j)}(\delta _{nj,1}\phi _{l_2,m_2}^{\prime \prime }(R)+{\displaystyle \frac{4\pi }{(2l_2+1)}}{\displaystyle _0^R}dr^{}D_{l_2}(R,r^{})\delta \rho _{l_2,m_2}^{v(nj)}(r^{}))+`$
$`+\delta \rho _{l_1,m_1}^{v(j)}(R)\left(\delta _{nj,1}\phi _{l_2,m_2}^{}(R)+{\displaystyle \frac{4\pi }{(2l_2+1)}}{\displaystyle _0^R}𝑑r^{}G_{l_2}(R,r^{})\delta \rho _{l_2,m_2}^{v(nj)}(r^{})\right)+`$
$`+\delta _{nj,1}k^2\delta \rho _{l_1,m_1}^{s(j)}\phi _{l_2,m_2}(R)\}`$
$`{\displaystyle \frac{e}{mR^2}}{\displaystyle \underset{j=1}{\overset{n1}{}}}{\displaystyle \frac{1}{nj}}{\displaystyle \underset{\stackrel{l_1,m_1}{l_2,m_2}}{}}I_2(l,m|l_1,m_1,|l_2,m_2)\times `$
$`\times \delta \rho _{l_1,m_1}^{s(j)}\left(\delta _{nj,1}\phi _{l_2,m_2}(R)+{\displaystyle \frac{4\pi }{(2l_2+1)}}{\displaystyle _0^R}𝑑r^{}r^2B_{l_2}(R,r^{})\delta \rho _{l_2,m_2}^{v(nj)}(r^{})\right)`$ (47)
The set of non-linear equations (46) and (47) must be solved iteratively starting from $`n=1`$. It is clear from the form of these equations that the index $`n`$ corresponds to the order of perturbation theory on the external field $`\phi (𝐫)`$. It is also seen that the chain of equations for the volume plasmon oscillations is self-sufficient, while the solutions for the surface oscillations also depend on those for volume, which is physically clear, because the electric field caused by volume oscillations of electron density must influence surface oscillations.

From (46) and (47), one derives the following solutions for $`n=1`$:

$$\delta \rho _{l,m}^{v(1)}(r)=\frac{\omega _p^2}{4\pi }\frac{k^2\phi _{l,m}(r)}{\omega ^2\omega _p^2}$$
(48)
$$\delta \rho _{l,m}^{s(1)}=\frac{\omega _p^2}{4\pi }\frac{1}{\omega ^2\omega _l^2}\left\{\phi _{l,m}^{}(R)\frac{4\pi (l+1)}{(2l+1)}_0^R𝑑r^{}\frac{r^{l+2}}{R^{l+2}}\delta \rho _{lm}^{v(1)}(r^{})\right\}$$
(49)
These equations show that, if the external field (30) is characterized by a certain angular momentum $`l_o`$ and its projection $`m_o`$, i.e. $`\phi _{l,m}(r)=\delta _{l,l_o}\delta _{m,m_o}\phi _{l_o,m_o}(r)`$, then $`\delta \rho _{l,m}^{s(1)}=\delta _{l,l_o}\delta _{m,m_o}\delta \rho _{l_o,m_o}^{s(1)}`$ and $`\delta \rho _{lm}^{v(1)}(r)=\delta _{l,l_o}\delta _{m,m_o}\delta \rho _{l_o,m_o}^{v(1)}(r)`$. Assuming these dependencies, the solutions of (46) and (47) for $`n=2`$ read as:

$`\delta \rho _{l,m}^{v(2)}(r)={\displaystyle \frac{e}{m((2\omega )^2\omega _p^2)}}I_1(l,m|l_o,m_o,|l_o,m_o)\times `$
$`\times \{\delta \rho _{l_o,m_o}^{v(1)}(r)(\phi _{l_o,m_o}^{}(r)+{\displaystyle \frac{4\pi }{(2l_2+1)}}{\displaystyle _0^R}dr^{}G_{l_o}(r,r^{})\delta \rho _{l_o,m_o}^{v(1)}(r^{}))`$
$`k^2\delta \rho _{l_o,m_o}^{v(1)}(r)\phi _{l_0,m_o}(r)4\pi \delta \rho _{l_o,m_o}^{v(1)}(r)\delta \rho _{l_o,m_o}^{v(1)}(r)\}`$
$`{\displaystyle \frac{e}{mr^2((2\omega )^2\omega _p^2)}}I_2(l,m|l_o,m_o,|l_o,m_o)\times `$
$`\times \delta \rho _{l_o,m_o}^{v(1)}(r)\left(\phi _{l_o,m_o}(r)+{\displaystyle \frac{4\pi }{(2l_o+1)}}{\displaystyle _0^R}𝑑r^{}r^2B_{l_o}(r,r^{})\delta \rho _{l_o,m_o}^{v(1)}(r^{})\right)`$ (50)
$`\delta \rho _{l,m}^{s(2)}={\displaystyle \frac{4\pi eN(l+1)}{m(2l+1)VR^{l+2}((2\omega )^2\omega _l^2)}}{\displaystyle _0^R}𝑑r^{}r^{l+2}\delta \rho _{l_o,m_o}^{v(1)}(r^{})+`$
$`+{\displaystyle \frac{e}{m((2\omega )^2\omega _l^2)}}I_1(l,m|l_o,m_o,|l_o,m_o)\times `$
$`\times \{\delta \rho _{l_o,m_o}^{s(1)}(\phi _{l_o,m_o}^{\prime \prime }(R)+{\displaystyle \frac{4\pi }{(2l_o+1)}}{\displaystyle _0^R}dr^{}D_{l_o}(R,r^{})\delta \rho _{l_o,m_o}^{v(1)}(r^{}))+`$
$`+\delta \rho _{l_o,m_o}^{v(1)}(R)\left(\phi _{l_o,m_o}^{}(R)+{\displaystyle \frac{4\pi }{(2l_o+1)}}{\displaystyle _0^R}𝑑r^{}G_{l_o}(R,r^{})\delta \rho _{l_o,m_o}^{v(1)}(r^{})\right)+`$
$`+k^2\delta \rho _{l_o,m_o}^{s(1)}\phi _{l_o,m_o}(R)\}`$
$`{\displaystyle \frac{e}{mR^2((2\omega )^2\omega _l^2)}}I_2(l,m|l_o,m_o,|l_o,m_o)\times `$
$`\times \delta \rho _{l_o,m_o}^{s(1)}\left(\phi _{l_o,m_o}(R)+{\displaystyle \frac{4\pi }{(2l_o+1)}}{\displaystyle _0^R}𝑑r^{}r^2B_{l_o}(R,r^{})\delta \rho _{l_o,m_o}^{v(1)}(r^{})\right)`$ (51)
By performing similar transformations, one can find the solutions $`\delta \rho _{lm}^{v(n)}(r)`$ and $`\delta \rho _{lm}^{s(n)}`$ for arbitrarily large $`n`$, although the formulae become more and more tedious the larger $`n`$ becomes. These formulae demonstrate that, in the higher orders of perturbation theory, plasmon resonances with angular momenta larger than the angular momentum of the external field can be excited. Indeed, the selection rules for the integrals $`I_1(l,m|l_o,m_o,|l_o,m_o)`$ and $`I_2(l,m|l_o,m_o,|l_o,m_o)`$ (see B) show that the angular momentum in $`\delta \rho _{lm}^{v(2)}(r)`$ and $`\delta \rho _{lm}^{s(2)}`$ can be twice as large as $`l_o`$. Equations (46\- 51) also demonstrate that the plasmon resonances in $`\delta \rho _{lm}^{v(n)}(r)`$ and $`\delta \rho _{lm}^{s(n)}`$ arise when $`\omega =\omega _p/n`$ and $`\omega =\omega _l/n`$ respectively. These equations indicate a significant shift of the plasmon resonance profiles towards lower frequencies in the highest orders of perturbation theory.

These results have a simple physical explanation. Absorption of several quanta of the external field (photons) by the cluster leads to the excitation of non-dipole plasmon oscillations of the electron density.

### 3.5 Fast electron-cluster collisions

Equations (46\- 51) can be used for the analysis of the balance between the surface and volume plasmon oscillations in the cluster. We demonstrate this for the example of fast electron scattering on a metal cluster. In this case, the external field of the projectile electron can be characterized by the Fourier component of the Coulomb potential

$$\phi (r)=\frac{4\pi }{q^2}e^{i\mathrm{𝐪𝐫}},$$
(52)
where $`𝐪=𝐩𝐩^{}`$ is the transferred momentum of the scattered electron.

The partial expansion of this potential reads as:

$$\phi (r)=4\pi \underset{l=0}{\overset{\mathrm{}}{}}\underset{m=}{\overset{m=l}{}}i^l\phi _{l,m}(r)Y_{l,m}^{}(𝐧_q),$$
(53)
where the partial component of the potential $`\phi _{l,m}(r)`$ is equal to:

$$\phi _{l,m}(r)=\frac{4\pi }{q^2}j_l(qr)Y_{l,m}(𝐧_r),$$
(54)
and $`j_l(qr)`$ is the spherical Bessel function (for definition see e.g.).

The form and properties of $`\phi _{l,m}(r)`$ \- (54) are exactly the same as assumed in (19) and (30). Therefore, from (48) and (49), one can immediately derive:

$$\delta \rho _{l,m}^{v(1)}(r)=\frac{\omega _pj_l(qr)}{\omega ^2\omega _p^2}$$
(55)
and

$$\delta \rho _{l,m}^{s(1)}=\frac{\omega _p}{q^2(\omega ^2\omega _l^2)}\left\{j_l^{}(qR)\frac{q^2\omega _p(l+1)}{(\omega ^2\omega _p^2)(2l+1)}_0^R𝑑r^{}\frac{r^{l+2}}{R^{l+2}}j_l(qr^{})\right\}$$
(56)
In the case of inelastic electron scattering, $`\omega `$ has the meaning of the transferred energy in the collision, $`\omega =\mathrm{\Delta }\epsilon =\epsilon \epsilon ^{}`$. Calculating the integral in (56) with the use of the well known properties of spherical Bessel functions (see e.g. ), one derives

$`\delta \rho _{l,m}^{s(1)}`$ $`=`$ $`{\displaystyle \frac{(2l+1)\omega _l^2}{\omega ^2\omega _l^2}}{\displaystyle \frac{j_l(qR)}{q^2R}}{\displaystyle \frac{\omega _p^2}{\omega ^2\omega _p^2}}j_{l+1}(qR)`$ (57)
Expresssions (55) and (57) coincide with those calculated in  in the plasmon resonance approximation by purely electrodynamic means, as the response of a dielectric sphere, having dielectric permeability $`ϵ=1\theta (Rr)\omega ^2/\omega _p^2`$.

From (55) and (57), one can easily elaborate the electron inelastic scattering cross section in the plasmon resonance approximation, using the method described in :

$`{\displaystyle \frac{d^2\sigma }{d\epsilon ^{}d\mathrm{\Omega }}}`$ $`=`$ $`{\displaystyle \frac{4p^{}R}{\pi pq^4}}{\displaystyle \underset{l}{}}(2l+1)^2j_l^2(qR){\displaystyle \frac{\omega _l^2\mathrm{\Delta }\epsilon \mathrm{\Gamma }_l^s}{(\mathrm{\Delta }\epsilon ^2\omega _{l}^{}{}_{}{}^{2})^2+\mathrm{\Delta }\epsilon ^2\mathrm{\Gamma }_l^{s2}}}+`$
$`+`$ $`{\displaystyle \frac{2p^{}R^3}{\pi pq^2}}{\displaystyle \underset{l}{}}(2l+1){\displaystyle \frac{\omega _p^2\mathrm{\Delta }\epsilon \mathrm{\Gamma }_l^v}{(\mathrm{\Delta }\epsilon ^2\omega _{p}^{}{}_{}{}^{2})^2+\mathrm{\Delta }\epsilon ^2\mathrm{\Gamma }_l^{v2}}}\times `$
$`\times `$ $`\left(j_l^2(qR)j_{l+1}(qR)j_{l1}(qR){\displaystyle \frac{2}{qR}}j_{l+1}(qR)j_l(qR)\right)`$
This cross section is totally determined by collective electron excitations in the cluster. The first and the second terms in (3.5) describe the contributions of the surface and the volume plasmon excitations respectively. In (3.5), we have also indroduced widths, $`\mathrm{\Gamma }_l^s`$ and $`\mathrm{\Gamma }_l^v`$, of the surface and volume plasmon resonances. They originate from the Landau damping mechanism of the plasmon excitations. For their determination we refer to the recent paper .

### 3.6 Multiphoton absorption

Next, we apply equations (46) and (47) to the description of the multiphoton absorption process. In this paper, we focus our consideration on the analysis of plasmon excitations. If surface or volume plasmon resonances are excited by photons, i.e. $`\omega \omega _p`$, then it is easy to check that the following condition is fulfilled $`\omega R/c\omega _pR/c1`$, where $`c`$ is the velocity of light. This condition implies the validity of the dipole approximation.

In the dipole approximation, one can neglect the momentum of the photon and put $`k=0`$. In this case, equations (46) and (47) are simplified dramaticaly. Indeed, from (46), one derives

$$\delta \rho _{l,m}^{v(n)}(r)=0$$
(59)
This result also simplifies equation (47) significantly for $`\delta \rho _{lm}^{s(n)}`$. After some trivial transformations it reduces to:

$`\left((\omega n)^2\omega _l^2\right)\delta \rho _{l,m}^{s(n)}={\displaystyle \frac{Ne^2}{mV}}\delta _{n,1}\phi _{lm}^{}(R)+`$
$`+{\displaystyle \frac{e}{m}}{\displaystyle \underset{j=1}{\overset{n1}{}}}{\displaystyle \frac{\delta _{nj,1}}{nj}}{\displaystyle \underset{\stackrel{l_1,m_1}{l_2,m_2}}{}}I_1(l,m|l_1,m_1,|l_2,m_2)\delta \rho _{l_1,m_1}^{s(j)}\phi _{l_2,m_2}^{\prime \prime }(R)`$
$`{\displaystyle \frac{e}{mR^2}}{\displaystyle \underset{j=1}{\overset{n1}{}}}{\displaystyle \frac{\delta _{nj,1}}{nj}}{\displaystyle \underset{\stackrel{l_1,m_1}{l_2,m_2}}{}}I_2(l,m|l_1,m_1,|l_2,m_2)\delta \rho _{l_1,m_1}^{s(j)}\phi _{l_2,m_2}(R)`$ (60)
The partial component of the linearely polarized dipole photon field is equal to

$$\phi _{l,m}(r)=\sqrt{\frac{4\pi }{3}}Er\delta _{l,1}\delta _{m,0}$$
(61)
Here $`E=\sqrt{2\pi \mathrm{}\omega /V_o}`$ is the strength of the photon’s electric field and $`V_o`$ is the normalization volume of the photon mode. Substituting (61) into (47), one derives

$`\left((\omega n)^2\omega _l^2\right)\delta \rho _{l,m}^{s(n)}=\sqrt{{\displaystyle \frac{4\pi }{3}}}{\displaystyle \frac{Ne^2E}{mV}}\delta _{n,1}\delta _{l,1}\delta _{m,0}+`$
(62)
$`+\sqrt{{\displaystyle \frac{4\pi }{3}}}{\displaystyle \frac{eE}{mR}}{\displaystyle \underset{l_1,m_1}{}}I_2(l,m|l_1,m_1,|1,0)\delta \rho _{l_1,m_1}^{s(n1)}`$ (63)
This equation should be solved iteratively starting from $`n=1`$. For $`n=1`$, the single non-trivial solution, $`\delta \rho _{1,0}^{s(1)}`$, reads as

$$\delta \rho _{1,0}^{s(1)}=\sqrt{\frac{4\pi }{3}}\frac{Ne^2E}{mV(\omega ^2\omega _1^2)}$$
(64)
Then, for $`n=2`$, the solution of (63) is of the form

$`\delta \rho _{l,m}^{s(2)}=\sqrt{{\displaystyle \frac{4\pi }{3}}}{\displaystyle \frac{eE\delta \rho _{1,0}^{s(1)}}{mR((2\omega )^2\omega _l^2)}}I_2(l,m|1,0|1,0)`$ (65)
The selection rules for $`I_2(l,m|1,0|1,0)`$ (see B) show that this integral does not vanish, when $`l=0`$ and $`m=0`$ or $`l=2`$ and $`m=0`$. Therefore, for $`n=3`$ from (63), one derives

$`\delta \rho _{l,m}^{s(3)}=\sqrt{{\displaystyle \frac{4\pi }{3}}}{\displaystyle \frac{eE}{mR((3\omega )^2\omega _l^2)}}`$
$`\{I_2(l,m|0,0|1,0)\delta \rho _{0,0}^{s(2)}+I_2(l,m|2,0|1,0)\delta \rho _{2,0}^{s(2)}\}`$ (66)
In (66), only the second term in brackets gives a non-zero contribution, since $`I_2(l,m|0,0|1,0)=0`$ (see B). Substituting $`\delta \rho _{0,0}^{s(2)}`$ and $`\delta \rho _{2,0}^{s(2)}`$ from (65) and $`\delta \rho _{1,0}^{s(1)}`$ from (64) in (66) and using the explicit expressions for the angular integrals given in B, one obtains

$`\delta \rho _{0,0}^{s(2)}={\displaystyle \frac{\pi ^{1/2}}{3m^2RV}}{\displaystyle \frac{Ne^3E^2}{\omega ^2(\omega ^2\omega _1^2)}}`$
$`\delta \rho _{2,0}^{s(2)}={\displaystyle \frac{4\pi ^{1/2}}{3\sqrt{5}m^2RV}}{\displaystyle \frac{Ne^3E^2}{(\omega ^2\omega _1^2)((2\omega )^2\omega _2^2)}}`$
$`\delta \rho _{1,0}^{s(3)}=\left({\displaystyle \frac{64\pi }{3}}\right)^{1/2}{\displaystyle \frac{16+3\sqrt{5}}{75m^3R^2V}}{\displaystyle \frac{Ne^4E^3}{(\omega ^2\omega _1^2)((2\omega )^2\omega _2^2)((3\omega )^2\omega _1^2)}}`$
$`\delta \rho _{3,0}^{s(3)}=\left({\displaystyle \frac{4\pi }{7}}\right)^{1/2}{\displaystyle \frac{12(2+\sqrt{5})}{75m^3R^2V}}{\displaystyle \frac{Ne^4E^3}{(\omega ^2\omega _1^2)((2\omega )^2\omega _2^2)((3\omega )^2\omega _3^2)}}`$ (67)
## 4 Induced multipole moments in the cluster

Let us now calculate the multipole moments of the cluster induced by an external radiation field on the basis of the model developed in the previous section, and analyse their plasmon resonance structure.

The induced multipole moment of the cluster is equal to

$$Q_{l,m}=\sqrt{\frac{4\pi }{2l+1}}𝑑Vr^lY_{l,m}(𝐧_r)\delta \rho (𝐫)$$
(68)
where the variation of electron density $`\delta \rho (𝐫)`$ is determined in (29) and (42). Subsituting (29) and (42) in (68) and putting $`\delta \rho _{l,m}^{v(n)}(r)=0`$ for any $`n`$ in the dipole approximation as follows from (59), one derives

$$Q_{l,m}^{(n)}=\sqrt{\frac{4\pi }{2l+1}}R^{l+2}\delta \rho _{l,m}^{s(n)}$$
(69)
Substituting in (69) $`\delta \rho _{1,0}^{s(n)}`$ from (64), one obtains the expression for the dipole moment of the cluster induced in the single-photon absorption process

$$D^{(1)}(\omega )=Q_{1,0}^{(1)}=\frac{Ne^2E}{m(\omega ^2\omega _1^2+\mathrm{i}\omega \mathrm{\Gamma }_1)}$$
(70)
The explicit expressions for the partial electron density variations $`\delta \rho _{l,m}^{s(n)}`$ entering (69) for $`n=2`$ (two photon case) and $`n=3`$ (three photon case) are given in (67). Subsituting the partial electron density variations $`\delta \rho _{0,0}^{s(2)}`$ and $`\delta \rho _{2,0}^{s(2)}`$ from (67) into (69), one derives the expression for the monopole and quadrupole moment of the cluster induced in the two photon regime

$`M^{(2)}(\omega )=Q_{0,0}^{(2)}={\displaystyle \frac{1}{2m^2R^2}}{\displaystyle \frac{Ne^3E^2}{\omega ^2(\omega ^2\omega _1^2+\mathrm{i}\omega \mathrm{\Gamma }_1)}}`$
$`Q^{(2)}(\omega )=Q_{2,0}^{(2)}={\displaystyle \frac{2}{5}}{\displaystyle \frac{Ne^3E^2}{m^2(\omega ^2\omega _1^2+\mathrm{i}\omega \mathrm{\Gamma }_1)((2\omega )^2\omega _2^2+i2\omega \mathrm{\Gamma }_2)}}`$ (71)
Here, we have introduced the plasmon resonance widths $`\mathrm{\Gamma }_1`$ and $`\mathrm{\Gamma }_2`$ which take into account Landau damping of the dipole and quadrupole surface plasmon resonances. They must be determined separately, e.g. by an ab initio approach (see ).

By absorbing three photons one can induce dipole and octupole moments in the cluster. Substituting $`\delta \rho _{1,0}^{s(3)}`$ from (67) into (69), one derives the expression for the induced dipole moment

$`D^{(3)}(\omega )`$ $`=`$ $`Q_{1,0}^{(3)}={\displaystyle \frac{4(16+3\sqrt{5})}{75m^3R^2}}\times `$ (72)
$`\times `$ $`{\displaystyle \frac{Ne^4E^3}{(\omega ^2\omega _1^2+\mathrm{i}\omega \mathrm{\Gamma }_1)((2\omega )^2\omega _2^2+i2\omega \mathrm{\Gamma }_2)((3\omega )^2\omega _1^2+\mathrm{i3}\omega \mathrm{\Gamma }_1)}}`$
The expression for the octupole moment induced by 3 photons following from (67) and (69) reads as

$`O^{(3)}(\omega )`$ $`=`$ $`Q_{3,0}^{(3)}={\displaystyle \frac{12(2+\sqrt{5})}{175m^3}}\times `$ (73)
$`\times `$ $`{\displaystyle \frac{Ne^4E^3}{(\omega ^2\omega _1^2+\mathrm{i}\omega \mathrm{\Gamma }_1)((2\omega )^2\omega _2^2+i2\omega \mathrm{\Gamma }_2)((3\omega )^2\omega _3^2+\mathrm{i3}\omega \mathrm{\Gamma }_3)}}`$
Here, we have also introduced the octupole plasmon resonance width $`\mathrm{\Gamma }_3`$.

Expressions (70-73) demonstrate that the multipole moments induced in the cluster during multiphoton absorption processes possess a prominent plasmon resonance structure. The nature of these resonances is the same as occurs in the mutiphoton absorption cross sections discussed in section 2.

The connection between $`D^{(1)}(\omega )`$ from (70) and the cross section $`\sigma _1`$ found in (6) is straightforward

$$\sigma _1=\frac{4\pi \omega }{cE}\mathrm{𝐼𝑚}D^{(1)}(\omega )$$
(74)
In the multiphoton regime, the connection between the induced multipole moments of the cluster and the multiphoton absorption cross section becomes more complex, which is apparent from the the classical nature of the expressions (71-73) and the explicit dependence of the multiphoton absorption cross sections on Planck’s constant. The discussion of this interesting relationship is however beyond the scope of the present paper.

## 5 Conclusion

In this paper, we have developed a formalism which allows one to calculate the cross section for multiphoton absorption in the plasmon resonance approximation. We have demonstrated that plasmon excitations with angular momenta larger than 1 substantially alter the profile for multiphoton absorption as compared to the single-photon case.

Our model is formulated in terms of a charge density distribution function $`\rho (𝐫)`$ for the cluster. This means that, in principle, one can study the response for different charge density profiles including deformed ones. Our model is a semi-classical one, which neglects the granularity of charge in the system. This is consistent with the principles underlying the jellium picture. It is appropriate for metallic clusters and, to a lesser extent, for fullerenes.

In the classical formulation of our model, we have used Euler’s equation for hydrodynamic flow, together with the equation of continuity. We have demonstrated that the results following from our model are consistent with direct estimates of the matrix elements for the multiphoton absorption process. The theoretical formalism we have developed is not confined in its application to photons. It can also be used to describe any kind of higher order plasmon excitation processes, for example those which arise by multiple scattering of electrons within a cluster.

## 6 Acknowledgments

The authors acknowledge support from the Royal Society of London, INTAS, the DAAD and the Alexander von Humboldt Foundation.

## Appendix A  Matrix elements of plasmon resonance transitions

In this appendix, we evaluate the matrix elements of plasmon resonance transitions in the plasmon resonance approximation by the use of the sum rule.

For a stationary (i.e. explicitly independent of time) operator $`\widehat{F}`$ of an observable physical quantity, characterising a system of particles with the Hamiltonian $`\widehat{H}`$, one can formulate the following sum rule (see e.g. )

$$\underset{n}{}\omega _{n0}\left|n|\widehat{F}|0\right|^2=\frac{1}{2}0|[\widehat{F},[\widehat{H},\widehat{F}]]|0$$
(75)
Here $`\omega _{n0}=\epsilon _n\epsilon _0`$, the summation is performed over all excited states of the system and $`[\widehat{H},\widehat{F}]`$ denotes the commutator of the operators $`\widehat{H}`$ and $`\widehat{F}`$: $`[\widehat{H},\widehat{F}]=\widehat{H}\widehat{F}\widehat{F}\widehat{H}`$.

Applying the sum rule (75) to operator $`\widehat{F}`$ defined as

$$\widehat{F}=\underset{k}{}F(𝐫_k),$$
(76)
where $`F(𝐫_k)`$ is a function of the coordinates of the $`kth`$ electron and the summation in (76) is performed over all particles in the system, one derives

$$\underset{n}{}\omega _{n0}\left|n|\widehat{F}|0\right|^2=\frac{\mathrm{}^2}{2m}𝑑𝐫\left|F(𝐫)\right|^2\rho (𝐫)$$
(77)
Here $`\rho (𝐫)`$ is the ground state charge density distribution in the system. Applying the general rule (77) to the function

$$F(𝐫)=\sqrt{\frac{4\pi }{2l+1}}r^lY_{lm}(𝐧)$$
(78)
and performing the integration in (77) with the density distribution (41), one derives

$$\underset{n}{}\omega _{n0}\left|n|\widehat{Q}_{lm}|0\right|^2=\frac{\mathrm{}^2}{2}\omega _l^2R^{2l+1}$$
(79)
Here, we have used the following definition of the multipole moments operator

$$\widehat{Q}_{lm}=\sqrt{\frac{4\pi }{2l+1}}\underset{k=1}{\overset{N}{}}e_kr_k^lY_{lm}(𝐧_𝐤).$$
(80)
The plasmon resonance frequencies in (79) are defined according to (2).

Using (79), one can easily evaluate the matrix elements of plasmon resonance transitions in the plasmon resonance approximation. Indeed, assuming that plasmon excitations dominate in the sum over $`n`$ in (79), one derives

$$Q_{lm}=n|\widehat{Q}_{lm}|0=\sqrt{\frac{\mathrm{}\omega _lR^{2l+1}}{2}}=eR^{l1}\sqrt{\frac{N\mathrm{}}{2m\omega _l}}\sqrt{\frac{3l}{2l+1}}$$
(81)
Equation (81) gives the matrix elements of plasmon resonance transitions for an arbitraty large angular momentum $`l`$. The correctness of the result (81) can be independently varified by performing calculations of the multipole dynamic polarizabilty of the cluster in the plasmon resonance approximation. Indeed, using (81), one derives

$$\alpha _l(\omega )=2\underset{n}{}\frac{\omega _{n0}|Q_{l0}|^2}{\omega _{n0}^2\omega ^2i\omega \mathrm{\Gamma }_l}\frac{R^{2l+1}\omega _l^2}{\omega _l^2(\mathrm{}\omega )^2i\mathrm{}^2\omega \mathrm{\Gamma }_l}$$
(82)
which is the known expression for the dynamic multipole plasmon polarizability in the plasmon resonance approximation (see e.g. ).

For the dipole plasmon resonance transition, one derives from (81)

$$Q_{10}=d_{10}=ez_{10}=e\sqrt{\frac{N\mathrm{}}{2m\omega _1}}$$
(83)
which is consistent with the dipole sum rule (4).

Equation (81) allows one to evaluate the matrix elements of electronic transitions between various plasmon resonance states. To demonstrate this, let us rewrite (81) in the form

$$Q_{lm}=\sqrt{\frac{1}{2l+1}}𝑑rr^{l+2}\rho _{l0}(r)$$
(84)
Here, we have introduced the radial transition density $`\rho _{l0}(r)`$ between the ground state and the excited plasmon resonance state with angular momentum $`l`$ and used the relationship $`I_1(l_n,m_n|l,m|0,0)=\delta _{l_n,l}\delta _{m_n,m}/\sqrt{4\pi }`$, when calculating angular intergrals in (84).

The radial transition density $`\rho _{l0}(r)`$ is localised in the vicinity of the cluster surface. Qualitatively, this is clear because $`\rho _{l0}(r)`$ describes the plasmon excitation. Quantitatively, this was proved by ab initio computations of the transition densities in the $`Na_{40}`$ and $`Na_{92}`$ clusters within the jellium model in . Therefore, to a reasonable accuracy, one can approximate $`\rho _{l0}(r)`$ by the delta function

$$\rho _{l0}(r)=\rho _{l0}\delta (rR)$$
(85)
Substituting (85) into (84) and comparing the result of the calculation with (81), one can determine the value $`\rho _{l0}`$ entering (85). The result of this calculation reads as

$$\rho _{l0}=\frac{e}{R^3}\sqrt{\frac{3l\mathrm{}N}{2m\omega _l}}$$
(86)
Let us now evaluate the matrix element for the dipole transition between plasmon resonance modes, which reads as

$$z_{l_2l_1}=\sqrt{\frac{4\pi }{3}}𝑑𝐫Y_{1,0}(𝐧)r\rho _{l_2l_1}(𝐫)$$
(87)
Here, $`\rho _{l_2l_1}(𝐫)`$ is the electron transition density between the dipole and the quadrupole plasmon modes. This transition density can be evaluated via the transition densities $`\rho _{l_2,l_1}(r)`$, $`\rho _{l_20}(r)`$ and the ground state electron density of the cluster $`\rho _{00}=eN/V`$ as follows

$$\rho _{l_2l_1}(𝐫)=\frac{\rho _{l_20}(r)\rho _{l_10}(r)}{\rho _{00}}Y_{l_1,0}(𝐧)Y_{l_2,0}(𝐧)$$
(88)
Substituting (88) into (87) and performing simple transformations, one derives

$$z_{l_2l_1}=e\sqrt{\frac{4\pi l_1l_2}{3}}\left(\frac{l_1(2l_2+1)}{l_2(2l_1+1)}\right)^{1/4}I_1(l_2,m_2|1,0|l_1,m_1)\frac{2\pi \mathrm{}}{m\omega _1}\delta (0)$$
(89)
This equation has an uncertainty, which originates from the fact that we have assumed zero width for the domain in the vicinity of the cluster surface in which plasmon excitations take place. By introducing a finite width $`\mathrm{\Delta }R`$ for this domain and using one of the standard representations of the $`\delta `$-function  to resolve the uncertainty, $`\delta (0)2/\pi \mathrm{\Delta }R`$, one finaly derives

$$z_{l_2l_1}=8e\sqrt{\frac{\pi l_1l_2}{3}}\left(\frac{l_1(2l_2+1)}{l_2(2l_1+1)}\right)^{1/4}I_1(l_2,m_2|1,0|l_1,m_1)\frac{\mathrm{}}{m\omega _1\mathrm{\Delta }R}$$
(90)
The explicit expression for the angular intergral $`I_1(l_2,m_2|1,0|l_1,m_1)`$ is given in B. In the case of the transition between the dipole and quadrupole plasmon resonance states this integral is equal to $`I_1(2,0|1,0|1,0)=1/\sqrt{5\pi }`$. Substituting this value into (90) and performing simple algebraic transformations, one arrives at the expression for the matrix element describing the transition between the dipole and the quadrupole plasmon resonance modes.

$$z_{12}=e\frac{8}{3}\left(\frac{6}{5}\right)^{1/4}\frac{\mathrm{}}{m\omega _1\mathrm{\Delta }R}$$
(91)
## Appendix B  Integrals $`I_1(l,m|l_1,m_1,|l_2,m_2)`$ and $`I_2(l,m|l_1,m_1,|l_2,m_2)`$

The angular integral,

$$I_1(l,m|l_1,m_1|l_2,m_2)=d\mathrm{\Omega }_{𝐧_r}Y_{l,m}^{}(𝐧_r)Y_{l_1,m_1}(𝐧_r)Y_{l_2,m_2}(𝐧_r),$$
(92)
is well known and can be found in many textbooks (see e.g. ). It is equal to

$`I_1(l,m|l_1,m_1|l_2,m_2)`$ $`=`$ $`(1)^mi^{l_1+l_2l}\sqrt{{\displaystyle \frac{(2l+1)(2l_1+1)(2l_2+1)}{4\pi }}}\times `$ (93)
$`\times `$ $`\left(\begin{array}{ccc}l& l_1& l_2\\ m& m_1& m_2\end{array}\right)\left(\begin{array}{ccc}l& l_1& l_2\\ 0& 0& 0\end{array}\right)`$
Here, the integral is expressed via 3j-symbols (for definition see e.g. ).

The angular integral,

$`I_2(l,m|l_1,m_1|l_2,m_2)`$ $`=`$ $`\sqrt{l_1(l_1+1)l_2(l_2+1)}{\displaystyle }d\mathrm{\Omega }_{𝐧_r}Y_{l,m}^{}(𝐧_r)\times `$ (94)
$`\times 𝐘_{l_1,m_1}^{(1)}(𝐧_r)𝐘_{l_2,m_2}^{(1)}(𝐧_r),`$
can be expressed via the sum of products of 3j-symbols and 6j-symbols (for definition see e.g. ) after performing the following transformations.

Using the standard relationships for spherical vector harmonics, written in  on page 210, one derives

$`I_2(l,m|l_1,m_1|l_2,m_2)=\sqrt{l_1(l_1+1)l_2(l_2+1)}{\displaystyle }d\mathrm{\Omega }_{𝐧_r}Y_{l,m}^{}(𝐧_r)\times `$
$`\times (\sqrt{{\displaystyle \frac{l_1+1}{2l_1+1}}}𝐘_{l_1,m_1}^{l_11}(𝐧_r)+\sqrt{{\displaystyle \frac{l_1}{2l_1+1}}}𝐘_{l_1,m_1}^{l_1+1}(𝐧_r))\times `$
$`\times \left(\sqrt{{\displaystyle \frac{l_2+1}{2l_2+1}}}𝐘_{l_2,m_2}^{l_21}(𝐧_r)+\sqrt{{\displaystyle \frac{l_2}{2l_2+1}}}𝐘_{l_2,m_1}^{l_2+1}(𝐧_r)\right)`$ (95)
The integrations arising in (95) can be performed and expressed via the sum of products of 3j-symbols and 6j-symbols, using the standard formulae given in  on pages 222 and 236. The result of the calculations of these integrals reads as

$`I_2(l,m|l_1,m_1|l_2,m_2)=`$
$`=(1)^{l_1+l_2+l+m+1}\sqrt{l_1(l_1+1)l_2(l_2+1)(2l+1)}\left(\begin{array}{ccc}l_1& l_2& l\\ m_1& m_2& m\end{array}\right)\times `$
$`\times \{\sqrt{{\displaystyle \frac{(l_1+1)(l_2+1)(2l_11)(2l_21)}{4\pi }}}\left(\begin{array}{ccc}l_11& l_21& l\\ 0& 0& 0\end{array}\right)\left\{\begin{array}{ccc}l_11& l_21& l\\ l_2& l_1& 1\end{array}\right\}+`$ (98)
$`+\sqrt{{\displaystyle \frac{l_1(l_2+1)(2l_1+3)(2l_21)}{4\pi }}}\left(\begin{array}{ccc}l_1+1& l_21& l\\ 0& 0& 0\end{array}\right)\left\{\begin{array}{ccc}l_1+1& l_21& l\\ l_2& l_1& 1\end{array}\right\}+`$ (101)
$`+\sqrt{{\displaystyle \frac{(l_1+1)l_2(2l_11)(2l_2+3)}{4\pi }}}\left(\begin{array}{ccc}l_11& l_2+1& l\\ 0& 0& 0\end{array}\right)\left\{\begin{array}{ccc}l_11& l_2+1& l\\ l_1& l_2& 1\end{array}\right\}+`$ (104)
$`+\sqrt{{\displaystyle \frac{l_1l_2(2l_1+3)(2l_2+3)}{4\pi }}}\left(\begin{array}{ccc}l_1+1& l_2+1& l\\ 0& 0& 0\end{array}\right)\left\{\begin{array}{ccc}l_1+1& l_2+1& l\\ l_2& l_1& 1\end{array}\right\}\}`$ (107)
For the particular cases of interest, one derives from (107)

$`I_2(1,0|0,0|1,0)`$ $`=`$ $`0`$
$`I_2(0,0|1,0|1,0)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{\pi }}}`$
$`I_2(2,0|1,0|1,0)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{5\pi }}}`$
$`I_2(1,0|2,0|1,0)`$ $`=`$ $`{\displaystyle \frac{6}{10\sqrt{\pi }}}\left(1+{\displaystyle \frac{16}{3\sqrt{5}}}\right)`$
$`I_2(3,0|2,0|1,0)`$ $`=`$ $`{\displaystyle \frac{18}{5\sqrt{105\pi }}}\left(1+{\displaystyle \frac{\sqrt{5}}{2}}\right)`$ (108)
## References


format: text

text:
# 1 Introduction

## 1 Introduction

Reconstruction and study of charmonium mesons in $`B`$ decays is a crucial component of the measurement of time-dependent CP-violating asymmetries.

The analyses described in the following paper are based on a sample of 20.7 fb-1 collected with BABAR at the $`\mathrm{{\rm Y}}(4S)`$resonance with an additional 2.6 fb-1 collected below the $`B\overline{B}`$threshold. A determination of the $`B`$ meson branching fractions depends upon an accurate measurement of the number of $`B`$ mesons in the data sample. The number of $`B\overline{B}`$events is determined by comparing the rate of multi-hadron events in data collected both on and off resonance. The continuum contribution to the on-resonance sample is estimated by rescaling the number of off-resonance hadronic events by the ratio of the number of observed $`\mu ^+\mu ^{}`$ events in the two samples. This procedure yields a total of $`22.72\pm 0.36`$ million $`B\overline{B}`$events.

## 2 Inclusive decays of $`B`$ to states containing Charmonium

$`J/\psi `$ candidates are selected by requiring two identified leptons of opposite charge. Electrons are identified based on the ratio $`E/p`$ of the energy deposited in the calorimeter to the measured momentum from tracking information, the shape of the calorimetric cluster and the ionization in the tracking detectors. Muons are identified by requiring a minimum ionizing signal in the calorimeter. In addition, the shape and penetration of the distribution of hits in the instrumented flux return are used. The number of $`J/\psi `$ events is determined by fitting the invariant mass distribution to a probability density function obtained from a simulation including contribution from both final state radiation and bremsstrahlung. The fit yields $`15739\pm 177`$ $`J/\psi `$ $`e^+e^{}`$ and $`13683\pm 154`$ $`J/\psi `$ $`\mu ^+\mu ^{}`$ signal events (Figure 1).

The $`\psi (2S)`$candidates are reconstructed in the $`\psi (2S)`$$`l^+l^{}`$ and $`\psi (2S)`$$`J/\psi \pi ^+\pi ^{}`$ decays. In the latter case the signal yield is extracted by a fit to the mass difference between the $`\psi (2S)`$and the $`J/\psi `$ reconstructed candidates. We find $`552\pm 50`$ decays to $`e^{}e^{}`$, $`437\pm 44`$ decays to $`\mu ^+\mu ^{}`$, $`405\pm 37`$ decays to $`J/\psi (e^+e^{})\pi ^+\pi ^{}`$ and $`400\pm 34`$ decays to $`J/\psi (\mu ^+\mu ^{})\pi ^+\pi ^{}`$.

The $`\chi _{c1}`$and $`\chi _{c2}`$candidates are selected by combining a reconstructed $`J/\psi `$ with a photon. The signal yield is determined by fitting the mass difference between the $`\chi _c`$ and the $`J/\psi `$ candidates. The shape of the signal is extracted from Monte Carlo, with the mass difference between the $`\chi _{c1}`$and $`\chi _{c2}`$peaks fixed to the PDG value . The fit gives $`476\pm 71`$ $`\chi _{c1}`$and $`86\pm 59`$ $`\chi _{c2}`$candidates for the $`J/\psi e^+e^{}`$ decay and $`545\pm 60`$ $`\chi _{c1}`$and $`104\pm 56`$ $`\chi _{c2}`$candidates for the $`J/\psi \mu ^+\mu ^{}`$ decay.

Values for the branching fractions are extracted from the yields. A $`90\%`$ CL limit is set on the $`B`$ decay to $`\chi _{c2}`$. The results for the inclusive production of $`J/\psi `$, $`\psi (2S)`$and $`\chi _c`$ in $`B`$ decays are summarized in Table 1, where the first uncertainty is statistical error and the second is systematics.

## 3 Exclusive decays of $`B`$ to Charmonium

The reconstruction of exclusive decay modes containing charmonium presents a very low background in most of the channels. For this reason the lepton identification criteria are loosened for one of the two $`J/\psi `$ decay products. As in the inclusive analysis, $`\psi (2S)`$candidates are reconstructed by their decay to $`e^{}e^{}`$, $`\mu ^+\mu ^{}`$ and $`J/\psi \pi ^+\pi ^{}`$. The $`\chi _c`$ candidates are selected through their decay to $`J/\psi \gamma `$.

The charmonium states are selected in a window around their expected mass for the decays to leptons. In the decays to $`J/\psi `$ states the mass difference distribution between the charmonium candidate and the reconstructed $`J/\psi `$ is used instead.

Selected candidates are then paired with a $`K^+`$, $`K_s^0`$(either $`\pi ^+\pi ^{}`$or $`\pi ^0\pi ^0`$), $`K^+`$ (either $`K^+`$$`\pi ^0`$ or $`K_s^0`$$`\pi ^+`$), $`K^0`$(either $`K^+\pi ^{}`$ or $`K_s^0`$$`\pi ^0`$), $`\pi ^0`$ or $`K_L`$ to form a $`B`$ candidate. The two most significant observables used to identify the signal are the difference $`\mathrm{\Delta }E`$in the center-of-mass frame between the reconstructed $`B`$ energy and half the nominally available energy, $`\sqrt{s}/2`$, and the energy-substituted mass, $`m_{ES}=\sqrt{s/4P_B^2}`$, where $`P_B^{}`$ is the center-of-mass momentum of the $`B`$ candidate. A sample of these distributions is given in Figure 2 for $`J/\psi K_S`$ events.

In the case of multiple candidates per event, only the candidate with the smallest $`|\mathrm{\Delta }E|`$ is selected. For all modes except $`BJ/\psi K_L`$ and $`BJ/\psi K^{}`$, the number of signal events is determined from the observed number of events in the $`(\mathrm{\Delta }E,m_{ES})`$ region after background subtraction. In addition to the usual combinatorial component, the background distribution shows an excess of events in the signal region. The combinatorial contribution is estimated by using an ARGUS function in the fit to the $`m_{ES}`$ distribution. The peaking background component is obtained from simulation of inclusive $`B`$ decays to charmonium, after removing the signal events.

The signal yields for the $`K^0`$ and $`K^+`$ modes are determined simultaneously from a likelihood fit, which is needed to account for the cross-feed between the $`K^{}`$ decay channels.

A different technique is used for the $`BJ/\psi K_L`$ decay mode. In this case only the $`K_L`$ direction is measured with information from the calorimeter and the instrumented flux return. Given this direction and the reconstructed charmonium candidate, the $`K_L`$ energy is extracted by using the $`B`$ mass as a constraint. To eliminate cross-feed from other decay modes, a veto has been introduced for events which have been selected already in other exclusive modes. This procedure yields a purity of about $`50\%`$. Due to this method, the $`m_{ES}`$ distribution cannot be used to determine the signal yield. The $`\mathrm{\Delta }E_{J/\psi K_L^0}`$ distribution is used in a log-likelihood fit. The shapes of the signal and inclusive charmonium background components are taken from Monte Carlo simulations. The shape of the non-charmonium background component is taken from an ARGUS fit to the $`\mathrm{\Delta }E_{K_L^0}`$ distribution for events in the $`J/\psi `$ mass sideband. After the background subtraction, this channel gives a signal yield of $`183\pm 14`$ events (Figure 3, left plot).

The information on the signal yields and the measured branching fractions for all presented exclusive modes is summarized in Table 2, while Figure 4 shows a comparison of the new preliminary BABAR results to the current PDG values.

In addition to the analyses described above, a very preliminary study has been recently performed by BABAR, on exclusive reconstruction of $`B`$ decays to $`\eta _cK`$ modes. The $`\eta _c`$ has been studied in the $`K_SK^+\pi ^{}`$, $`K^+K^{}\pi ^0`$ and $`K^+K^{}K^+K^{}`$ ($`\varphi \varphi `$ plus non-resonant) decay modes. The energy-substituted mass distribution for the decay $`\eta _cK_SK^+\pi ^{}`$ is presented as an example in Figure 3, right plot. No value has been extracted yet for the branching fractions in these decay modes.

## 4 Summary

Using 22.7 million $`B\overline{B}`$events recorded by the BABAR detector, the inclusive branching fractions for the production of $`J/\psi `$, $`\psi (2S)`$and $`\chi _c`$ are presented. Combining the charmonium state with either a $`K^\pm `$, $`K^0`$, $`K^\pm `$, $`K^0`$or $`\pi ^0`$, $`B`$ decays are reconstructed exclusively and their branching fractions are determined.