Surface photoeﬀect and optical absorption: Maxwell equations in tensor and nonlocal isotropic approximation

The surface photoeﬀect is studied for a p polarized continuous laser with an angular frequency ω in the long wave length (LWL) domain, i.e. λ > 100 ˚A. In this LWL domain the laser ﬁelds are usually taken to be constant in space, an approximation called LWL approximation. However, perpendicular to a gas-solid interface, the potential, the electron density and the laser ﬁelds vary on a sub-nanometric scale. Here, these varying laser ﬁelds are solutions of the combined Amp`ere-Maxwell and material equations. The material functions are taken to be either scalar or tensor. They are a function of the unperturbed electron density of the material system calculated from the Schr¨odinger equation. For a planar structureless gas-solid interface and a p polarized laser, the Amp`ere-Maxwell equation in real space reduces to a set of coupled ﬁrst and second order ordinary diﬀerential equations in one dimension. For an Al(001) surface calculation using model jellium and DFT/LDA potentials, the total absorbed power density of the laser and the number of electrons escaping from the surface, display monopole and multipole surface plasmon resonances and a minimum near the bulk plasmon resonance found by other calculations and in experiments. [DOI: 10.1380/ejssnt.2009.249]


I. INTRODUCTION
The surface photoeffect at a planar gas-solid interface in the long wave length domain (LWL), i.e. λ 100Å or, equivalently, for energies ω 124 eV, has been widely discussed in the past (see e.g. the review papers by Kliewer [1] and Feibelman [2]). One concluded that, at the gas-solid interface and in the bulk, the LWL approximation, corresponding to spatially constant laser fields, does not apply.
The first reliable experimental study using the photon angle and energy resolved photoyield (AERPY) with the electron ejected normal to the surface of the clean electrochemically polished Al(001) surface was performed by Levinson, Plummer and Feibelman [3,4]. The initial states were the Fermi level and the surface state. This spectrum shows a maximum at about 0.8ω p =12.5 eV, interpretted later to correspond to the multiple surface plasmon resonance, and a minimum at the bulk plasmon frequency ω p =15.3 eV. Using the model of Feibelman [2,5,6] based on a jellium potential, these authors reproduced the photoelectron yield of the Al(001) surface. More recently, similar photoelectron spectrum experiments were performed by Barman et al. [7] on a clean Al(111) surface.
The surface photoelectric effect can be identified in the AERPY and electron energy loss (EELS; see e.g. Tsuei et al. [8]) experimental spectra either by the presence of an asymmetric line shape of the direct excitation or of a multipole surface plasmon resonance. In AERPY the photoelectron yield Y is | f | E · ∇|i + f | ∇· E|i | 2 , where E is the electric field of the laser. The first term is called standard, direct or transverse and the second term surface, indirect or longitudinal. If the electric field is taken independent of the spatial coordinate, only the direct term survives. Therefore the second, surface term is responsible for the asymmetric line shapes and multipole surface plasmon resonances. The surface or indirect excitations have been used to explain the experimental spectra of simple metals like Na, K, Cs (see ref. [8]) and Al (see above) and transition metals like Cu(001) [9], Ag(111) [10] and Ag(001) [11,12].
In the standard theory of optical reflection at a planar interface the imposed continuity of the displacement D implies a discontinuity of the associated electric field E of the laser. To eliminate this discontinuity, one either calculates the corrections to the ordinary reflection/refraction optics formulas (see e.g. Feibelman [2], Apell [13] and Kempa and Schaich [14]) or directly solves the Ampère-Maxwell equation. To solve this equation, methods can be divided into macroscopic and microscopic, and classical, semi-classical and quantum. The widely used ones are the semi-classical infinite barrier method (see the review of Kliewer [1]), the infinite barrier method (see e.g. Gerhardts and Kempa [15]), the hydrodynamic method (see e.g. Kempa et al. [16]), the method of Mukhopadhyay and Lundqvist [17] and the quantum microscopic method of Feibelman (see the review ref. [2]). In the majority of these methods an Al(001) gas-solid system was studied, where the tensor and non local material functions were obtained using a jellium solid and the Ampère-Maxwell equation was solved in momentum space. The cited authors use tensor and non local material functions in their models and think that (see e.g. Kliewer [1] and Feibelman [2]) simpler scalar and local material functions will disagree with the experimental results.
Here the Ampère-Maxwell equation in real space is written in terms of the vector potential A in a semiclassical (classical laser field, quantum electron density) approximation. The standard temporal gauge of solid state physics is used, i.e. the scalar potential is set to zero. The Ampère-Maxwell equation is combined with the material equations where the susceptibility of the actual gas-solid system is taken to be a tensor and a nonlocal isotropic function. This susceptibility is a function of i) the analytic susceptibilities of a jellium bulk dependent on the laser frequency and on the photon wave vector q and ii) the unperturbed electron density of the material system calculated by solving the Schrödinger equation. The associated model is called the vector potential from electron density (VPED). As the gas-solid interface is seen by a laser in the LWL domain as a structureless plane, the laser field parallel to the surface r can be taken to be a plane wave. Then the Ampère-Maxwell equation for a p polarized laser reduces to a set of coupled first and second order ordinary differential equations (CODE) in one dimension perpendicular to the surface. The VPED-CODE model is solved directly through fully numerical algorithms. As a first application of the model, we study a classical example namely the absorbed power density of the laser by a Al(001) surface. We will compare the absorbed power density calculated from VPED using either scalar and local or tensor and nonlocal material functions. These calculations can not be compared directly to the AERPY or EELS spectra but we will be able to identify in the photoabsorption spectrum the monopole and multipole surface plasmon resonances and the minimum associated with the bulk plasmon resonance.

II. MAXWELL AND MATERIAL EQUATIONS
For a linearly polarized laser incident on a planar gassolid interface one defines a plane of polarization or of incidence (POI; here xOz) by the laser incident direction q and the direction z normal to the surface. For the present gas-solid system there are two coordinate systems of interest: the photon coordinate system defined by the propagation direction of the light q with the photon electric and magnetic fields perpendicular to it and the surface coordinate system defined by the above POI and a plane parallel to the surface. In this paper we use the temporal gauge in the surface coordinate system, i.e. the scalar potential set to zero (see a very pedagogical discussion of the gauges used in physics by Heras [18]). We restrict the discussion to a p linearly polarized laser. Then the electric field E of the laser in the surface coordinate system has two projections: x transverse or standard and z longitudinal or surface. If the contribution of the longitudinal component to an observable is not negligible then the surface photoeffect appears.
We solve the semi classical Ampère-Maxwell equation involving the current density written in SI or in the atomic unit system, i.e.
where we restrict the discussion to non magnetic materials with the permeability µ = µ 0 . We combine the above equation with the following material equations written in tensor form where ε 0 is the vacuum dielectric function,1 is the identity tensor and "·" is the dot product. In the above equatioñ χ b ,χ c andσ are the relative bound and conduction susceptibility and the conductivity tensors. Note that these tensors are complex and include the absorption of the photon by the material system. To obtain the material functions needed in the Ampère-Maxwell equation (Eq. (1)) for the classical electromagnetic fields our vector potential from the electron density (VPED) model performs the following steps: 1. The bulk conduction susceptibility tensorχ c is written in momentum space as a product of two factors: the scalar plasma frequency ω p and a tensor factor F function of the variables q and ω. The plasma frequency ω p is itself written as a constant times a "scale factor", in fact the electron density ρ s of the material system over a unit volumẽ In the photon coordinate system,F is a second rank 2×2 diagonal tensor with transverse or standard F t and longitudinal or surface F l components. Following Linhard [19], Kliewer and Fuchs [1,20] and Mermin [21] the transverse component F t in the jellium approximation reads where q = q 2 x + q 2 z is the module of the wave number of the photon, u = q/(2 k F ) and v = (ω + iΓ)/(q v F ). The Fermi velocity, the Fermi wave number and the relaxation time of the electrons are designated by v F , k F and 1/Γ. In the same approximation, the longitudinal F l component reads 2. For the gas-solid system the bulk electron density ρ s is replaced by ρ(z, ω) ρ s . The relative electron density ρ(z, ω) normalized to ρ s is the result of a calculation where the bulk, the surface and the vacuum regions are taken into account. Due to the computational constraints that require zero susceptibility at the boundaries of the calculation grid a thick slab including two vacuum regions is used in our modeling. The relative electron density ρ(z, ω) is proportional to the square of the wave functions Here ψ i ( r) and φ i (z) are the i-th unperturbed wave functions of the material system in three-and onedimension spaces. The one-dimensional discrete states φ i (z) are obtained in a box of dimension L after solving the Schrödinger equation using the discrete variable representation (DVR) method developed by Light et al. [24]. Here we use the formulation of the DVR proposed by Colbert and Miller [25]. In eq. (8), the factor k 2 F − k 2 i corresponds to the number of plane wave functions or the density of states of the conduction electrons parallel to the surface. The unperturbed electron density of an Al(001) surface is presented in Fig. 1. To calculate this density two e-matter model potentials were used: a jellium potential in the Jennings et al. [23] parametrization with the Wigner-Seitz radius of r s = 2 a.u. and a model potential fitted by Chulkov et al. [22] to a density functional theory (DFT) calculation using a local density approximation (LDA) for the exchange-correlation functional. Compared to the jellium, the DFT/LDA potential presents a shifting of the surface edge and periodic oscillations due to the presence of an effective realistic valence electronic structure. One can shift the jellium edge to fit to DFT/LDA model potential rise but in the present work we use the standard non shifted jellium potential;

The Fourier transform of the electron density ρ(z, ω)
in momentum space gives ρ(q z , ω). Note that our conjugate variables are z and q z and from the laws of optics q x is conserved over all space; 4. The conduction electrons tensor susceptibilitỹ χ c ( q, ω) is multiplied first by the momentum space electron density ρ(q z , ω), rotated from the photon to surface coordinate system using matrix T (see e.g. Kliewer and Fuchs [20]) and finally converted to the real space by an inverse Fourier transform where q = q 2 x + q 2 z and in the above expression we used Eq. (3) with ρ s replaced by ρ(q z , ω) ρ s . Because q, q z and q x /q are variables of the above Fourier integral, the matrix T is also an implicit variable. In other words, following Kliewer and Fuchs [20], one defines the angle ϕ = tan −1 (q x /q z ) an implicit variable in the Fourier integral that does not equal θ i the fixed incident angle of the photon; 5. The scalar bulk bound electrons susceptibility χ b s is obtained from the analysis of the experimental data (for Aluminium we used Ehrenreich and Phillip [26]). Then this susceptibility is simply multiplied by the electron density in real space. In the calculation of the Linhard tensor bulk susceptibility (Eqs. (4), (5), (7) and (6)) for the conduction electrons we again used the data of Ehrenreich and Phillip (numerical values are given in results section). Using Eq. (9), the final total electron susceptibility of the bound and conduction electrons reads Note thatχ b (z, q x , ω) is a diagonal tensor with identical elements ρ(z, ω) χ b s (ω) on the diagonal. The obtained total susceptibilityχ(z, q x , ω), which will be used in the Ampère-Maxwell equations below, is a non diagonal complex susceptibility tensor in the surface coordinate system. Note that for the conduction electron susceptibility, the multiplication in momentum space ρ(q z , ω)F ( q, ω) corresponds to a convolution in direct spacẽ As it is well known, in the solid this convolution corresponds to a nonlocal isotropic case and therefore we designate the material functions of the present model tensor and nonlocal isotropic material functions.
Our model can immediately be reduced to a scalar and local form simply by writing the tensor susceptibility as where the scalar bulk susceptibility χ s (ω) can be either calculated from the Drude model for the conduction and obtained from the experiment for bound electrons (for Aluminium see Ehrenreich and Philipp [26]) or taken from the tables of Palik [27,28]. Again the relative electron density ρ(z, ω) is calculated from Eq. (8).
Models including the gas-solid interface in the material functions through the use of the electron density of the real system have been used in the past. For example Bagchi et al. [29] have studied a W(001) surface using a scalar and local model including a z dependent electron density whereas Apell [30] and Georges [31] used a tensor model for a jellium surface with analytic dispersive susceptibilities (our Eqs. (3), (4) and (6)) but solved the Ampère-Maxwell equation in the momentum space.
Next in the Ampère-Maxwell equation the laser fields of the photon are replaced by the vector potential A. In the temporal gauge, where the scalar potential is set to zero, one defines the vector potential A of the laser through its connection to the electric field E and the magnetic induction B Then, replacing Eqs. (2) and (13) For a planar interface between two media and a continuous laser (notation cw) of frequency ω, the vector potential can be written as a product of two terms where again q x is the x projection of the photon wave vector q, a constant for planar interfaces. Introducing Eq. (15) in (14) and performing the tem-poral derivatives, one obtains in the SI unit system where the total susceptibilityχ of equation (10) will be used below. For p polarization of the laser, the above equation generates a system of two coupled equations corresponding respectively to the unit vectors j x and j z . This is obtained by writing the scalar and dot products in terms of their components, multiplying from the left by exp(−i q x x + i ω t) and integrating over r and time.
The resulting system of coupled first and second order ordinary differential equations (CODE) for a linearly polarized laser p in one dimension z perpendicular to the surface reads This system of equations contains diagonal and off diagonal components of the susceptibility tensor. It can be used with the tensor susceptibility functions (Eq. (10)) where χ zx = 0 or with the scalar susceptibility functions (Eq. (12)) where χ zx = 0. In both cases one has a system of coupled equations. For the numerical calculation, the system of CODE is converted to a system of three first order differential equations, the derivatives discretized and solved numerically using the box method (see e.g. Zwillinger [32]; p.532) with the second order Adams-Multon or Brugnano and Trigiante [33] propagators. Full details of the numerical algorithm will be given in a forthcoming publication. Briefly we use a thick slab where the laser is incident from one side of the slab and this incident flux is preserved in the calculated laser field. This permits us to obtain the reflexion coefficient p r and the impedance p Z (see Kliewer [1]) where χ is a scalar susceptibility.

III. POWER DENSITY, ELECTRON ESCAPING PROBABILITY AND REFLECTANCE
The most comprehensible discussion of the power absorption including dispersion in the classical hydrodynamic model framework has been given by Forstmann et al. [16,34,35]. Following this author, the standard power density per unit volume averaged over t, x and y can be written as a sum of three terms where ρ = ∇ · P and −β 2 = B ε 0 ω 2 p (B is given by Frostmann). Here P is the polarization and X and Y define the size of the averaged box. The term proportional to Γ is an absorption term whereas the term proportional to β 2 , appearing because of the tensor nature of the dielectric functions, corresponds to an oscillatory exchange between the electromagnetic field and the material system. Forstmann argued that the evaluation of the power density absorption has to be restricted to the term proportional to Γ where (χ) is the imaginary part of the susceptibility and the squares correspond to the product of complex conju-gates. Following Kempa and Forstmann [16], the number http://www.sssj.org/ejssnt (J-Stage: http://www.jstage.jst.go.jp/browse/ejssnt/) of electrons reaching the surface can be written using the mean power density per volume unit, Eqs. (21) a n s (θ i , ω,ζ(ω)) = a Π(θ i , ω, z) exp(−z/ζ(ω)) t,X, where we recall thatζ(ω) is the mean escape length of the electron excited by the laser taken, following Penn [36], as a function of the laser frequency. I 0 is the mean incident laser intensity per time unit. Following Born and Wolf [37] (their sec. 1.5.3) the reflectance R per surface unit or the reflected/incident energy flux is proportional to the normal component S z of the Poynting vector. In the vacuum region on the laser side of the slab the S z can be separated in two components where S + z is the incident and S − z the reflected Poynting vector projections perpendicular to the surface. In the photon coordinate system the Poynting vector is proportional to the square of the laser electric field whereas in the surface coordinate system one writes The reflectance is the ratio of these reflected S − z and incident S + z components of the Poynting vector In the next section we will present the results of the power density absorption a Π ∼ J · J * /I 0 function of the penetration depth z in the solid, of the spectrum of the integrated power density absorption a Π ∼ Γ J · J * and of the spectrum of the number of escaping electrons a n(θ i , ω,ζ). We will also present the reflectance spectrum R. All these observables will be calculated using the scalar and tensor susceptibilities.

IV. RESULTS
The semi-classical model VPED-CODE presented in this paper will now be applied to the power density absorption of a p polarized cw laser incident at 45 • on a Al(001) a metal surface. The system of coupled Ampère-Maxwell equations in one dimension z (Eqs. (17)) is solved in real space. The susceptibility χ(z, q x , ω) is either a tensor, Eq. (10), or a scalar, Eq. (12). The scalar susceptibility of the conduction electrons was calculated using the standard Drude theory. In the expressions Eqs. (3), (4) and (6), we used the parameters of bound and conduction electrons of Al as given by Kittel and Ehrenreich, Phillips and Segal [26], namely ω p = 12.7 eV, v F = 2.02 10 6 m/s and τ = 1/Γ = 5.13 10 −15 s. The electron density (Eq. (8)), entering in the expression of susceptibility, is calculated using a slab of the material system 800Å thick with 400Å of vacuum on each side of the slab. In the studied frequency domain of the laser the skin depth of Al(001) is smaller than 280Å. For our slab thickness of 800Å less than 1/10 of the initial laser intensity reaches the outer edge of the slab. The wave functions used to calculate the unperturbed electron density are solutions of the Schrödinger equation with either a jellium potential (parametrization of Jennings et al. [23]) or a more realistic potential fitted by Chulkov et al. [22] to a DFT/LDA calculation. Only this last potential explicitly  introduces the electron-nuclei interaction of the atoms, takes into account the electron-electron and electron-hole interactions and accommodates the surface and the image states. Figure 2 presents the total a Π (Eq. (21)) power density absorption function of the penetration coordinate z in the spectral region of the bulk plasmon resonance (15.3 eV experiment, 15.7 eV theory). The results in the first and the second columns are obtained respectively using the jellium and DFT/LDA potentials. The absorption in the bulk rises slightly with increasing laser frequency. The surface contribution to the photoabsorption has a minimum at the bulk plasma frequency. Therefore the minimum in the power density absorption spectra (Figs. 3 and 4 below) is due to the minimum in the surface power density absorption for both jellium and Chulkov DFT/LDA potentials. Figures 3 and 4 display the spectra of the integrated power density absorption a Π (Eqs. (21), (22), (23)) and of the number of escaping electrons a n (Eq. (24)) using the jellium and the DFT/LDA potentials, respectively. These spectra can be compared qualitatively with the AERPY and EELS experimental spectra discussed in the introduction. But a precise comparison is impossible because, for example, the AERPY experiment measures the photoelectron yield Y ∼ | f | A · ∇|i + f | ∇ · A|i | 2 , from precise initial |i and final |f states. In the power density absorption and electron escaping probability one excites electrons from all the accessible states at a given laser frequency. The unperturbed electron density was calculated using the DFT/LDA potential of Chulkov et al. [22]. Same explanation as in Fig. 3.
The main trends of the presented spectra are similar for the two potentials of the Al(001) material system using both scalar and tensor susceptibilities. For the power density absorption one has two minima, two maxima and one inflexion at about 18 eV. In the escaping probability spectrum the features are less pronounced. The feature around 10.5 eV is probably related to the monopole and multipole surface resonances (see e.g. Barman et al. [7] AERPY of Al(111) and Tsuei et al. [8] for EELS spectra); the minimum at around 15 eV is related to bulk plasmon resonance; our feature at 18 eV probably corresponds to the two features in the experimental spectrum of Al(111) discussed by Barman et al. [7]. These authors interpret them as originating from the single and collective electron excitations. At low energy our power density absorption and electron escaping probability spectra display high intensities. In AERPY experiment of Levinson and Plummer [3,4], to be compared in first approximation to the electron escaping probability, the spectrum starts at 9 eV then rises at the multipole surface plasma frequency (see Fig. 6 below). Our calculated electron escaping probability does not fall strongly enough below 9 eV. The power density absorption spectrum can not be directly compared with AERPY because it excites all the accessible states at a given photon frequency. This power density absorption depends on the susceptibilities which are large at low laser frequencies.
From Figs. 3 and 4 it appears that the partial contributions a Π x and a Π z (Eqs. (22) and (23))) or a n x and a n z show opposite slopes. The minimum at the bulk plasma frequency corresponds to the crossing between these two contributions. When using scalar susceptibility this crossing and minimum seems to be shifted to lower frequencies. Figure 5 presents the spectrum of the reflectance R (Eq. tentials. The curve labelled "optics" corresponds to a calculation of R as the square of the reflexion coefficient p r (Eq. (18)) obtained directly from the material constants function of the laser frequency only. For a given potential, jellium (upper graph) or DFT/LDA (lower graph), there is almost no difference between the scalar and tensor susceptibilities. The calculations using the jellium potential show a reflectance very close to the result using material constants. The calculations using the model DFT/LDA potential for the scalar reflectance R, are close to the "optics" result but display a minimum at the plasma frequency. The tensor reflectance R, is also close to the "optics" result and presents no minimum at the bulk plasma frequency. For the moment we have no explanation for this minimum in the scalar reflectance though it can be related to the parametrization of the DFT/LDA potential.

V. DISCUSSION AND CONCLUSION
This paper presents a semi-classical model for the photon absorption at the gas-solid interface valid in the LWL spectral domain, i.e. λ 100Å. The Ampère-Maxwell equation, written for the vector potential of a p polarized photon, is solved in the classical approximation. This equation uses scalar or tensor material functions (see Eqs. (3), (10) and (9)) obtained as a product of the bulk material constants and the gas-solid slab electron density. The gas-solid slab electron density is calculated in the quantum approximation from the Schrödinger equation. In this equation, a jellium or a DFT/LDA (Chulkov et al. [22]) potential was used. Our model, called vector potential from electron density (VPED), solves a system of two coupled ordinary differential equations (CODE).
As discussed in the introduction, many models have been developed to solve the Maxwell equation at a planar interface. Our VPED-CODE model generates a set of coupled differential equations of the type obtained in the model of Mukhopadhyay and Lundqvist [17]. However, our equations are not integro-differential and we solve them in real space. Furthermore, instead of eliminating the second derivative as Mukhopadhyay and Lundqvist did, we numerically solve a system of three coupled ordinary differential equations of the first order.
In the power density absorption function of the penetration coordinate z (Fig. 2) the main difference between the jellium and the DFT/LDA potential calculations is that the latter potential displays larger oscillations of the power density absorption in the bulk corresponding to the oscillations in the electron density (see Fig.  1). But these oscillations influence only slightly the integrated power density absorption and electron escaping probabilities (Figs. 3 and 4) discussed below.
One of the goals of the present study was to compare the calculations of power density absorption for different potentials. The spectra of the integrated power density absorption calculated using the jellium and the DFT/LDA potentials (Figs. 3 and 4) display a similar behavior. The small quantitative differences observed in our results suggest that one does not need sophisticated potentials to obtain acceptable results for power density absorption at least for a simple free electron metal like aluminium. Figure 5 of Kempa and Forstmann [16] presented the power density absorption for aluminium at different laser frequencies calculated using several models. Many of the discussed models reproduce the minimum near the bulk plasma frequency but do not display the surface plasmon maximum at 12.5 eV and the feature at high laser frequencies. Our Figs. 3 and 4 reproduce these features slightly shifted to lower energies for all the presented calculations.
The power density absorption displays similar features with the photoelectron AERPY (Levinson, Plummer and Feibelman [3,4] and Barman et al. [7]) and electron energy loss EELS (Tsuei et al. [8]) spectra. But, as mentioned above, the observables associated with AERPY and EELS differ from the power density absorption in two aspects: the yield corresponds to a different operator and the initial and final states are well defined in AERPY whereas in power density absorption all the accessible initial states are excited. For example the experiment by Levinson, Plummer and Feibelman [3,4] measures the photoelectron spectrum corresponding to the excitation only from the Fermi and the surface states. Figure 6 displays our recent calculations [38] using the VPED-CODE model with scalar material functions of the photoelectron spectrum of aluminium together with experimental results of Levinson and Plummer [4]. Our model reproduces rather well the experimental photoelectron spectrum except for the multipole plasmon resonance at 12.5 eV. Further investigations with tensor material and non local susceptibilities are needed to verify this preliminary result.
Another goal of the present study was to compare the results with scalar and tensor material functions. The results using scalar (first line of Figs. 3 and 4) and tensor (second line of these figures) material functions are also very similar thus leading us to believe that simple scalar material functions can be used in the calculations at least at an exploratory level.
The Al(001) system studied here is known to be a free electron metal and was used as a benchmark to test the theoretical models. We have obtained rather satisfactory results by introducing bound and conduction electron contributions to the susceptibility. The calculations of other metallic and non metallic surfaces is possible. For the scalar material functions the calculations are simpler and do not need any preliminary analysis in terms of bound and conduction electrons. These calculations can be done using the material constants of the tables of Palik [27,28]. For the tensorial material constants the distinction between bound and conduction electrons allows a study of other transition metals where the analysis in terms of bound and conduction electrons contributions to the susceptibility is known. Such examples are silver and copper (see Ehrenreich and Phillip [39]).