{'fraction_non_alphanumeric': 0.039263635114713556, 'fraction_numerical': 0.016623162522324497, 'mean_word_length': 4.5728066146072575, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Formalism for Multiphoton Plasmon Excitation in Jellium Clusters

We present a new formalism for the description of multiphoton plasmon excitation processes in jellium clusters. By using our method, we demonstrate that, in addition to dipole plasmon excitations, the multipole plasmons (quadrupole, octupole, etc) can be excited in a cluster by multiphoton absorption processes, which results in a significant difference between plasmon resonance profiles in the cross sections for multiphoton as compared to single-photon absorption. We calculate the cross sections for multiphoton absorption and analyse the balance between the surface and volume plasmon contributions to multipole plasmons.


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.
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 [20] The inelastic scattering of fast electrons on metal clusters in the range of transferred energies above the ionization threshold was considered in [21]. 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 [21].
The role of the polarization interaction and plasmon excitations in the process of electron attachment to metal clusters has also been examined both theoretically [22,23] and experimentally [24]. 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 Appendix A, we show how matrix elements for collective transitions can be calculated on the basis of the sum rules. In Appendix B, we present details of calculations of the angular integrals arising in the formalism outlined here.

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.

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 Here, e is the charge of electron, c is the velocity of light,h is Planck's constant, ω no = ε n − ε 0 is the electron excitation energy, ω 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. [10,17] and section 4 of this paper) Here V = 4πR 3 /3 is the cluster volume, where R = r o N 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 1 − 4eV , 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 [11,17]) 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 [4,5,6,7,10,11]), which means that the summation in the sum rule (see e.g. [25]) need be performed only over excited states in the vicinity of the Mie resonance. Now, assuming a Lorentzian distribution of width Γ 1 for the plasmon resonance states and replacing the delta function, δ(ω no −hω), in (1) by the profile (see e.g. [25]) one recovers the well-known expression for the single-photon absorption cross section (see e.g. [4,10]) The width Γ 1 is due to Landau damping. Its calculation for metal clusters is performed, for example, in [21].
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.

Two-photon absorption
In the dipole approximation, the two-photon absorption cross section is equal to 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 n m z nm z mō hω − ω mo + iδ Here, we have also introduced a dipole plasmon resonance width Γ 1 and used the sum rule (4) for the evaluation of the matrix elements for the dipole plasmon excitation |r 10 | 2 ≈ Nh/2mω 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ω ∼ ω 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, δ(ω no − 2hω) by a Lorentzian distribution of width Γ 2 (see e.g. [25]) By substituting (8) and (9) in (7), one derives 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 ω = ω 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 [26]. By estimating the radial component of the momentum of a single electron in a dipole and quadrupole plasmon oscillatory mode as p 1 ∼ mω 1 ∆R and p 2 ∼ mω 2 ∆R respectively, one derives Here, A is a dimensionless constant, of the order of one, ∆z 1 ∼h/p 1 and ∆z 2 ∼h/p 2 are the uncertainties relating to an electron in the dipole and quadrupole plasmon modes respectively, ∆R is the width of the domain in the vicinity of the cluster surface within which plasmon excitations take place. In Appendix A we prove the correctness of this estimate and demonstrate that matrix element z 21 is equal to By substituting (12) in (10), we obtain the final expression for the two-photon absorption cross section: where A = 8 We note that the cross section (13) depends explicitly on Planck's constanth, while the cross section (6) does not. The independence of (6) fromh is connected with the fact that plasmon oscillations are a purely classical effect, while the dependence of (13) onh 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 singleand 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 Γ 1 = ω 1 /4, Γ 2 = ω 2 /4, ∆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.

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  (6) and (13) and normalised per unit atom. The two-photon absorption profile is scaled by a factor 1/100. The scales are not identical for the two curves for reasons of definition of the cross sections in the single-and two-photon cases, but both are given in atomic units.
the n-photon absorption cross section has the following form: Here the amplitude M n is equal to 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 ω = ω 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 twophoton 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 singleparticle character. This follows, for instance, from the fact that with increasing l the wave length of the surface plasmon mode, 2πR/l, becomes smaller than the characteristic wave length of the delocalised electrons at the Fermi surface, 2πh/ √ 2m∆ε, where ∆ε 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.

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.

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 dv/dt with the total local electric field E acting on the density at the point (r, t). It has the following form: The electric field E 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: Here ϕ(r, 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 δρ(r, t).
We assume that the external potential ϕ(r, t) is the solution of the wave equation. Therefore, we can put where ϕ(r, t) satisfies the equation Here k = ω/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 where ρ o (r) is the electron density distribution in a free cluster without an external field and δρ(r, 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: ∂ρ(r, t) ∂t The simultaneous solution of equations (17), (20) and (21) with appropriate initial conditions and the initial distribution ρ o (r) allow one to determine the variation of electron density δρ(r, t) as well as its velocity v(r, t). We solve this problem by using a perturbative approach on the external field ϕ(r, t).

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 E ≪ mω 2 R/e is fulfilled. Substituting here the characteristic values ω ∼ 0.1, R ∼ 10, one derives E ≪ 0.1 in atomic units or E ≪ 5 · 10 8 B/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: By substituting these expansions into (17), (21) and performing simple transformations, one derives v n (r) = ie mnω δ n,1 ∇ϕ(r) iωnδρ n (r) + ∇ · ρ 0 (r)v n (r) + n−1 Here δ n,1 is the Kroneker symbol. One can exclude v n (r) from equation (25) by the substitution of (24) in (25). Performing this transformation with the simultaneous use of (19) and ∆|r − r ′ | −1 = −4πδ(r − r ′ ), one derives the following equation: 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 ρ o (r) and ∇ρ o (r) 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 ϕ(r).
Solving the set of equations (24) and (26) with ϕ(r) describing the dipole electronphoton 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 δρ 1 (r) as follows from the purely electrodynamical perturbative approach to the electron scattering problem [21].

Spherically symmetric case
Equations (24) and (26) are valid for an arbitrary shape of the initial distribution ρ o (r). 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 ρ o (r) = ρ 0 (r). This relationship allows one easily to exclude angular variables from equation (26). Using this relationship together with the partial expansion for δρ n (r) and ϕ(r), one derives Here, Y l,m (n 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 (n r ) and then integrated over dΩ nr . We also used the well known expansion (see e.g. [27]) where function B l (r, r ′ ) is defined as follows Here Θ(r ′ − r) is the step function. In equation (31) we have introduced the function G l (r, r ′ ) = r 2 ∂B l (r, r ′ )/∂r, which is of the the form G l (r, r ′ ) = l r l−1 r ′l−1 Θ(r ′ − r) − (l + 1) 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 δρ n−j (r) and ϕ(r) according to (29) and (30), which we denote as The gradient terms in the sum contain the integration of a spherical harmonic multiplied by the scalar product of the two vector spherical harmonics: This type of integral arises, when one expresses the gradients of the potential and the density according to (see e.g. [27]) Here Y l,m (n r ) are, respectively, the longitudinal and the transverse vector spherical harmonics, the definition of which one can find in [27]. We mention some properties of these vector hamonics: Y (−1) l,m (n r ) = n r Y l,m (n r ), Y (1) l,m (n r ) = ∇ Ω Y l,m (n r )/ l(l + 1) and Y (−1) l,m (n r ) · Y (1) l,m (n r ) = 0. The potential Φ n−j l 2 ,m 2 (r) in (38) is as follows 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 Appendix B.
Using the formulae written above, one can easily rewrite (31) in the following form 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).

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 Here N is the total number of delocalized electrons in the cluster volume V = 4πR 3 v /3. In this case it is natural to look for the solution of (40) by expressing it in the following form 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 v → R 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 When deriving (43), we have used the fact that Θ ′ (R − r) = −δ(R − r). Also we have introduced the function D l (r, r ′ ) = l 2 (l 2 − 1) r l 2 −2 r ′l 2 −1 Θ(r ′ − r) + (l 2 + 1)(l 2 + 2) 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 ω l is the same as in (2). The volume plasmon resonance frequency is equal to In equation (43), ω 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: 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 ϕ(r). 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: lo,mo (r). Assuming these dependencies, the solutions of (46) and (47) for n = 2 read as: By performing similar transformations, one can find the solutions δρ v(n) lm (r) and δρ s(n) lm 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 lm arise when ω = ω p /n and ω = ω 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.

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 where q = p − p ′ is the transferred momentum of the scattered electron.
The partial expansion of this potential reads as: where the partial component of the potential ϕ l,m (r) is equal to: and j l (qr) is the spherical Bessel function (for definition see e.g. [25]). The form and properties of ϕ l,m (r) -(54) are exactly the same as assumed in (19) and (30). Therefore, from (48) and (49), one can immediately derive:

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. ω ∼ ω p , then it is easy to check that the following condition is fulfilled ωR/c ∼ ω p R/c ≪ 1, 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)

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 where the variation of electron density δρ(r) is determined in (29) Substituting in (69) δρ s(n) 1,0 from (64), one obtains the expression for the dipole moment of the cluster induced in the single-photon absorption process The explicit expressions for the partial electron density variations δρ s(n) l,m entering (69) for n = 2 (two photon case) and n = 3 (three photon case) are given in (67). Subsituting the partial electron density variations δρ s(2) 0,0 and δρ s(2) 2,0 from (67) into (69), one derives the expression for the monopole and quadrupole moment of the cluster induced in the two photon regime 2,0 = − 2 5 Here, we have introduced the plasmon resonance widths Γ 1 and Γ 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 [21]). By absorbing three photons one can induce dipole and octupole moments in the cluster. Substituting δρ s (3) 1,0 from (67) into (69), one derives the expression for the induced dipole moment The expression for the octupole moment induced by 3 photons following from (67) and (69) reads as Here, we have also introduced the octupole plasmon resonance width Γ 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) (ω) from (70) and the cross section σ 1 found in (6) is straightforward 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.

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 singlephoton case.
Our model is formulated in terms of a charge density distribution function ρ(r) 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.
Applying the sum rule (A1) to operatorF defined aŝ where F (r k ) is a function of the coordinates of the k − th electron and the summation in (A2) is performed over all particles in the system, one derives n ω n0 n|F |0 2 =h 2 2m dr |∇F (r)| 2 ρ(r) (A3) Here ρ(r) is the ground state charge density distribution in the system. Applying the general rule (A3) to the function F (r) = 4π 2l + 1 r l Y lm (n) (A4) and performing the integration in (A3) with the density distribution (41), one derives n ω n0 n|Q lm |0 2 =h 2 2 ω 2 l R 2l+1 (A5) Here, we have used the following definition of the multipole moments operator Q lm = 4π 2l + 1 N k=1 e k r l k Y lm (n k ).
The plasmon resonance frequencies in (A5) are defined according to (2). Using (A5), 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 (A5), one derives Q lm = n|Q lm |0 = hω l R 2l+1 2 = eR l−1 Nh 2mω l 3l 2l + 1 (A7) Equation (A7) gives the matrix elements of plasmon resonance transitions for an arbitraty large angular momentum l. The correctness of the result (A7) can be independently varified by performing calculations of the multipole dynamic polarizabilty of the cluster in the plasmon resonance approximation. Indeed, using (A7), one derives α l (ω) = 2 n ω n0 |Q l0 | 2 ω 2 n0 − ω 2 − iωΓ l ≈ R 2l+1 ω 2 l ω 2 l − (hω) 2 − ih 2 ωΓ l (A8) which is the known expression for the dynamic multipole plasmon polarizability in the plasmon resonance approximation (see e.g. [10,17]). For the dipole plasmon resonance transition, one derives from (A7) Q 10 = d 10 = ez 10 = e Nh 2mω 1 (A9) which is consistent with the dipole sum rule (4). Equation (A7) allows one to evaluate the matrix elements of electronic transitions between various plasmon resonance states. To demonstrate this, let us rewrite (A7) in the form Here, we have introduced the radial transition density ρ 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) = δ ln,l δ mn,m / √ 4π, when calculating angular intergrals in (A10).
The radial transition density ρ l0 (r) is localised in the vicinity of the cluster surface. Qualitatively, this is clear because ρ 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 [19]. Therefore, to a reasonable accuracy, one can approximate ρ l0 (r) by the delta function ρ l0 (r) = ρ l0 δ(r − R) (A11) Substituting (A11) into (A10) and comparing the result of the calculation with (A7), one can determine the value ρ l0 entering (A11). The result of this calculation reads as Let us now evaluate the matrix element for the dipole transition between plasmon resonance modes, which reads as z l 2 l 1 = 4π 3 drY 1,0 (n)rρ l 2 l 1 (r) (A13) Here, ρ l 2 l 1 (r) is the electron transition density between the dipole and the quadrupole plasmon modes. This transition density can be evaluated via the transition densities ρ l 2 ,l 1 (r), ρ l 2 0 (r) and the ground state electron density of the cluster ρ 00 = eN/V as follows ρ l 2 l 1 (r) = ρ l 2 0 (r)ρ l 1 0 (r) ρ 00 Y l 1 ,0 (n)Y l 2 ,0 (n) (A14)