Ab initio theory of the Drude plasma frequency
AAb initio theory of the Drude plasma frequency
Bernardo S. Mendoza ∗ and W. Luis Moch´an Centro de Investigaciones en ´Optica, Le´on, Guanajuato, M´exico Instituto de Ciencias F´ısicas, Universidad Nacional Aut´onoma de M´exico,Av. Universidad s/n, Col. Chamilpa, 62210 Cuernavaca, Morelos, M´exico.
Abstract
We derive a theoretical expression to calculate the Drude plasma frequency ω p based on quantum me-chanical time dependent perturbation theory. We show that in general ω p should be replaced by a secondrank tensor, the Drude tensor D , which we express in terms of products of the velocity expectation valuesintegrated over the Fermi surface, and which is amiable to analytical and numerical evaluation. For theSommerfeld’s model of metals our expression yields the ubiquitous plasma frequency ω p = 4 πn e e /m e . TheDrude tensor takes into account the geometry of the unit cell and may be calculated from first principles forisotropic as well as anisotropic metallic systems. We present results for the noble metals, Cu, Ag, and Auwithout stress and subject to isotropic and uniaxial strains, and we compare the results to those availablefrom experiment. PACS numbers: a r X i v : . [ c ond - m a t . o t h e r] N ov . INTRODUCTION The Drude theory of metals, put forward back in 1900, includes among its successes a simpleexplanation that allows for approximate estimates of properties of metals whose full comprehensionrequired the development of the quantum theory of condensed matter. Drude’s theory of electricaland thermal conductivity was based on the kinetic theory gases applied to the conduction electronsof the system, while the the core electrons strongly bonded to the atomic nuclei are taken as aninert entity. The precise assumptions of the Drude model are clearly explained in textbooks ofsolid state theory like Ref. 3. Drude assumed that the electronic velocity distribution is given bythe Maxwell-Boltzmann distribution function, which leads to the wrong result for the specific heatof metals of 3 k B /
2, where k b is Boltzmann’s constant. It was Sommerfeld, almost 25 years afterthe publication of Drude’s results, that recognized the fact that Fermi-Dirac (FD) statistics wererequired for treating the electron correctly. As it turns out, Sommerfeld’s model is essentially theclassical electron gas model used by Drude, but using the FD distribution for the valence electrons.Albeit, one of the first success of Sommerfeld was to obtain the experimental linear behavior ofthe specific heat of metals.The Drude model yields the following well known result for the dielectric function (cid:15) ( ω ) of ametal, (cid:15) ( ω ) = 1 − ω p ω , (1)where ω p = 4 πn e e m e , (2)is known as the Drude plasma frequency, where n e , − e and m e are the number density, charge andfree mass of the conduction electrons, and ω is the frequency of the light that perturbs the electrons.Although Eq. (1) gives a precise description of the interaction of light with metals for photonswith energy below the threshold of electronic interband transitions, a fully quantum-mechanicalderivation would justify its successes.In this article, we derive a closed expression for ω p based on quantum mechanical time dependentperturbation theory. We show that in general ω p should be replaced by a second rank tensor, thatwe call the Drude tensor D , for which we derive a closed expression that is amiable to analyticevaluation in the case of free independent electrons, and to numerical evaluation for ab initio quantum mechanical calculations. The Drude tensor takes into account the geometry of the unitcell and may be calculated from first principles for isotropic as well as anisotropic metallic systems.2he article is organized as follows. In Sec. II, we present the theoretical formalism, showing themain expressions used to calculate the Drude tensor. In Sec. II A we derive ω p within the Drude-Sommerfeld model as a special case of our formalism, and in Sec. II B we present the procedure tonumerically evaluate the Drude tensor. In Sec. III we present results for the noble metals, Ag, Cuand Au without stress and subject to an isotropic and to an uniaxial strain, and we summarizeour findings in Sec. IV. II. ANALYTIC EXPRESSION FOR THE DRUDE PLASMA FREQUENCY
In order to derive the analytic expression for the Drude plasma frequency, we assume theelectrons may be described through an independent particle approximation, although we do allowfor many-body effects through an effective Hamiltonian that depends on all occupied states, as indensity functional theory. The electrons interact with an electromagnetic field which we assumeis a classical field. Thus we describe quantum mechanical matter interacting with classical fields.We neglect local field and excitonic effects. We write the one electron Hamiltonianˆ H ( t ) = ˆ H + ˆ H I ( t ) , (3)as the sum of an unperturbed effective time-independent Hamiltonian ˆ H that describes the in-teraction of an electron with the crystalline lattice and its effective interaction with the otherelectrons, and an interaction Hamiltonian ˆ H I ( t ) that describes the interaction of the electron witha time-dependent electromagnetic field. We describe the state of the system through the one elec-tron density operator ˆ ρ , with which we can calculate the expectation value of any single-particleobservable ˆ O as (cid:104) ˆ O(cid:105) = Tr( ˆ ρ ˆ O ) with Tr denoting the trace. Within the interaction picture thedensity operator evolves in time due to the interaction Hamiltonian according to i (cid:126) ddt ˆ ρ ( t ) = [ ˆ H I ( t ) , ˆ ρ ( t )] , (4)while the operators that correspond to all observables evolve through ˆ H according toˆ O ( t ) = ˆ U † ( t ) ˆ O s ( t ) ˆ U ( t ) , (5)where ˆ O s ( t ) is the same observable in the Schr¨odinger picture, given by ˆ O (0) for operators thatdo not depend explicitly on time, and ˆ U ( t ) = exp( − i ˆ H t/ (cid:126) ) (6)3s the non-perturbed unitary time-evolution operator. Assuming the field is turned on adiabatically,we may integrate (4) to yield ˆ ρ ( t ) = ˆ ρ + 1 i (cid:126) (cid:90) t −∞ dt (cid:48) [ ˆ H I ( t (cid:48) ) , ˆ ρ ( t (cid:48) )] , (7)where ˆ ρ is the unperturbed, time-independent equilibrium density matrix. We look for the stan-dard perturbation series solution, ˆ ρ ( t ) = ˆ ρ + ˆ ρ (1) ( t ) + ˆ ρ (2) ( t ) + . . . , where the superscript denotesthe order (power) with which each term depends on the perturbation ˆ H I ( t ). Since we are interestedonly in the linear response, we concentrate our attention on the 1-st order termˆ ρ (1) ( t ) = 1 i (cid:126) (cid:90) t −∞ dt (cid:48) [ ˆ H I ( t (cid:48) ) , ˆ ρ ] . (8)We will take our system as a solid described by a non-perturbed periodic Hamiltonian, whoseeigenfunctions are Bloch states, | m k (cid:105) , characterized by a band index m and a crystal momentum k . For ˆ H I ( t ) we take an interaction with an electromagnetic field with a wavelength much largerthan the crystal parameter. Thus, electronic transitions due to this interaction are vertical, i.e.,they conserve k .Within the dipole approximation, the interaction Hamiltonian in the Length Gauge is given by ˆ H I ( t ) = e ˆ r ( t ) · E ( t ) , (9)where ˆ r ( t ) = ˆ U † ( t ) ˆ r (0) U ( t ) is the position operator of the electron at time t and E ( t ) the timedependent perturbing classical electric field.From (8) we obtain the first order density matrix elements between Bloch states ρ (1) nm ( k ; t ) ≡ (cid:104) n k | ˆ ρ (1) ( t ) | m k (cid:105) = 1 i (cid:126) (cid:90) t −∞ dt (cid:48) (cid:104) n k | [ ˆ H I ( t (cid:48) ) , ˆ ρ ] | m k (cid:105) = − ie (cid:126) (cid:90) t −∞ dt (cid:48) e iω nm ( k ) t (cid:48) (cid:104) n k | [ ˆ r (0) , ˆ ρ ] | m k (cid:105) · E ( t (cid:48) ) , (10)where ω nm ( k ) ≡ ω n ( k ) − ω m ( k ) and E n ( k ) = (cid:126) ω n ( k )) are the unperturbed energy eigenvaluescorresponding to the stationary Schr¨odinger’s equation ˆ H | n k (cid:105) = E n ( k ) | n k (cid:105) . Notice that ˆ ρ hasmatrix elements (cid:104) n k | ˆ ρ | m k (cid:105) = δ nm f ( E n ( k )) , (11)with f the Fermi-Dirac distribution, which at the temperature T = 0 becomes f ( E n ( k )) = Θ( E F − E n ( k )) ≡ f n ( k ) , (12)4ith E F the Fermi energy of the system and Θ the unit step function. This defines the distributionsfunctions f n ( k ) in reciprocal space, one for each band.It is convenient to represent the position operator in coordinate space ˆ r (0) → r , when cal-culating its interband matrix elements and in reciprocal space ˆ r (0) → i ∇ k when calculating its intraband matrix elements, so that following Ref. 6, we can readily show that (cid:104) n k | [ˆ r , ˆ ρ ] | m k (cid:105) = f mn ( k ) r nm ( k ) + iδ nm ∇ k f n ( k ) , (13)where f nm ( k ) = f n ( k ) − f m ( k ) and r nm ( k ) = (cid:104) n k | r | m k (cid:105) . Notice that f nm ( k ) = 0 if n = m . Thewell-known commutator ˆ v = ˆ˙ r = 1 i (cid:126) [ ˆ r , ˆ H ] , (14)allows us to write the interband matrix element as r nm ( k ) = v nm ( k ) iω nm ( k ) ( n (cid:54) = m ) , (15)where ˆ v is the velocity operator related to the momentum operator by ˆ p = m e ˆ v .To obtain the optical linear response we look for the expectation value of the macroscopicpolarization density P , whose time derivative yields the current density, i.e., the expectation valueof the current operator. Thus, ∂∂t P = − e Ω Tr(ˆ ρ (1) ( t )ˆ v ( t )) , (16)with Ω the volume of the unit cell. Assuming an harmonic perturbation E ( t ) = E ( ω ) e − iωt , weobtain the tensorial relation P a ( ω ) = χ ab ( ω ) E b ( ω ) , (17)where χ ab ( ω ) is the linear susceptibility response tensor, where the superscripts a , b denote Carte-sian components, and we use Einstein convention for repeated indices. Using Eqs. (10)-(17), weobtain χ ab ( ω ) = ie (cid:126) ω (cid:88) mn (cid:90) BZ d k π v amn ( k ) (cid:18) f mn ( k ) r bnm ( k ) + iδ nm ∂ k b f n ( k ) ω nm ( k ) − ω (cid:19) = χ abe ( ω ) + χ abi ( ω ) , (18)where the sums over band indices and over the wavevector correspond to the trace and the latterwas replaced by the usual integral over the first Brillouin zone (BZ), we defined ∂ k b = ∂/∂k b to5implify our notation, and we identified the interband χ abe ( ω ) and intraband χ abi ( ω ) contributionsto the susceptibility as those that contain the first and second terms in the numerator of Eq. (18),those which contain the factor f nm ( k ) and the Kronecker’s delta δ nm respectively. Using Eq. (15)we write the interband contribution as χ abe ( ω ) = e (cid:126) ω (cid:88) mn (cid:90) BZ d k π f mn ( k ) ω nm ( k ) (cid:18) r amn ( k ) r bnm ( k ) ω nm ( k ) − ω (cid:19) . (19)The intraband contribution is χ abi ( ω ) = ie (cid:126) ω (cid:88) mn (cid:90) BZ d k π v amn ( k ) (cid:18) iδ nm ∂ k b f n ( k ) ω nm ( k ) − ω (cid:19) = e (cid:126) ω (cid:88) n (cid:90) BZ d k π v an ( k ) ∂ k b f n ( k ) , (20)where we denote v ann ( k ) as v an ( k ), the electron’s velocity for the n -th band at point k of the Brillouinzone. We recognize χ abe ( ω ) as the well known linear susceptibility amply discussed in the scientificliterature as well as in textbooks of physics; in the rest of the article, we only devote our attentionto χ abi ( ω ), that as we see now, leads to a closed and computationally amenable expression for theDrude frequency.We rewrite Eq. (20) as χ abi ( ω ) = − πω D ab , (21)where we define the Drude tensor as D ab = − πe (cid:126) (cid:88) n (cid:90) BZ d k π v an ( k ) ∂ k b f n ( k )= e π (cid:88) n (cid:90) BZ d k v an ( k ) v bn ( k ) δ ( E F − E n ( k )) . (22)Here, we used Eq. (12) for f n ( k ), employed the derivative of the Heaviside step function withrespect to its argument, d Θ( x ) /dx = δ ( x ), and identified the velocity of the electron in the n -thband from the dispersion of the corresponding energy ∇ k E n ( k ) = (cid:126) v n ( k ).For any two scalar valued functions F ( k ) and G ( k ) of the crystal momentum k , we can evaluate (cid:90) BZ d k F ( k ) δ ( G ( k )) = (cid:90) S d σ k F ( k ) |∇ k G ( k ) | , (23)which takes us from an integration over the first Brillouin zone with a Dirac’s delta function δ ( G ( k ))to a surface integral over that surface S within the first Brillouin zone for which G ( k ) = 0, with d σ k the differential element of area in reciprocal space. Notice that S may have zero, one or more6onnected components, so that we interpret the surface integral as implying a sum over all of them.Applying this result to Eq. (22) we obtain D ab = e π s (cid:88) n (cid:90) S n d σ k v an ( k ) v bn ( k ) |∇ k ( E F − E n ( k )) | = e π (cid:126) s (cid:88) n (cid:90) S n d σ k v an ( k ) v bn ( k ) v n ( k ) , (24)where S n is the contribution of the n -th band to the Fermi Surface (FS), defined by those Blochvectors k for which E n ( k ) = E F . This is the main result of our paper, as it allows the calculationof the Drude tensor, the intraband susceptibility and its contribution to the dielectric response forany metal from its electronic structure. We have included in Eq. (24) a spin-degeneracy factor s which allows us to ignore the spin in the state labels n for those systems that are spin-degenerate.Thus, we take s = 1 if the band index includes the spin, as would be the case when the spin-orbitcoupling is included in the calculation, and s = 2 when the spin is not explicitly included.For a semiconductor there are no bands that cross the Fermi Energy, and then the Drude tensoris identically zero, or in other words, there is no intraband contribution to the susceptibility. Fora metal, there is at least one band that crosses the E F , yielding a finite intraband contribution.The intraband dielectric tensor may be written written as (cid:15) abi = δ ab + 4 πχ abi = δ ab − D ab ω , (25)giving the characteristic 1 /ω divergence as ω → D ab has units of frequency squared and that it would agree with Eq. (1) for a Drude modelif we identify D ab = ω p δ ab with ω p given by Eq. (2). Thus, the interband contribution to thedielectric response given by Eqs. (24) and (25) can be seen a generalization of the Drude model.However, this dielectric response has an explicit tensorial character in contrast to that found instandard textbooks, which always describe a scalar dielectric response. The symmetry of thesystem determines which components of D ab are finite and how they relate among themselves.We may write the Drude tensor in a more familiar way by integrating by parts the first line ofEq. (22), to obtain D ab = 4 πe (cid:126) (cid:88) n (cid:90) BZ d k π f n ( k ) ∂ k b v an ( k ) , (26)and identify ∂ k b v an ( k ) = (cid:126) ( m ∗ n ( k )) − ab , (27)7rom Ref. 8, where m ∗ n ( k ) is the effective mass tensor. Then D ab = 4 πe (cid:88) n (cid:90) BZ d k π f n ( k )( m ∗ n ( k )) − ab = 4 πn e e (cid:104) ( m ∗ ) − ab (cid:105) , (28)where the sum may be restricted to the conduction bands. This result is similar to the usualDrude formula with n e the number density of conduction electrons but with the inverse electronicmass 1 /m e replaced by an average (cid:104) ( m ∗ ) − (cid:105) of the inverse effective mass tensor over the occupiedstates of the conduction bands. Although this result looks appealing, it hides the fact that onlythe electrons at the Fermi surface do contribute to the dynamics and the response of the system.Furthermore, the numerical integration over the full BZ in Eq. (28) would require a larger grid { k α } of k-points and would be more costly to evaluate than the numerical integral (24) only overthose k points close to the Fermi surface (see Sec. II B). Finally, evaluation of the inverse masstensor in Eq. (28) requires a larger computational cost and has a larger numerical uncertaintythan the evaluation of the velocity matrix elements in Eq. (24). A. Ideal Metal
Within the Sommerfeld theory of metals the conduction electrons have the dispersion relationof free electrons E = (cid:126) k / m e and fill up a Fermi sphere of radius k F , related to the numberdensity n e of conduction electrons through n e = k F / π . The speed of electrons at the Fermisurface is v F = (cid:126) k F /m e . Due to the spherical symmetry, D ab = 0 for a (cid:54) = b , and D xx = D yy = D zz = D aa / ≡ ω p , Then, from Eq. (24) ω p = 13 e π (cid:126) s (cid:90) k F d Ω v F = 4 π k F π e m e = 4 πn e e m e , (29)where we wrote d σ k = k F d Ω, with d Ω the differential element of solid angle. Thus ω p is the wellknown Drude value for the plasma frequency, however derived from Eq. (24) which comes from apurely quantum mechanical approach within time-dependent perturbation theory, instead of thetext-book derivation from the phenomenological model of Drude. Curiously, ω p is proportional tothe number density n e of conduction electrons, although only those electrons at the Fermi surfacecontribute to the integral in Eq. (24). B. Realistic metal
Eq. (24) is an elegant expression of the Drude tensor in terms of an integral over the Fermisurface of a tensor formed by products of velocity components divided by the speed v an v bn /v n .8owever, for actual applications to realistic models it is convenient to return to a volume integral.To that end, we realize that for any function g ( k ) we can write (cid:90) S n d σ k g ( k ) = (cid:90) BZ d k g ( k ) |∇ f n ( k ) | , (30)as can be immediately verified using Eqs. (12) and (23). Thus, we write Eq. (24) as D ab = e π (cid:126) s (cid:88) n (cid:90) BZ d k |∇ k f n ( k ) | v an ( k ) v bn ( k ) v n ( k ) . (31)To carry out the integration numerically, we generate a regular grid { k α } of points insidethe Irreducible Brillouin Zone that corresponds to the crystallographic group of the metal understudy, numbered by a discrete set of indices α . Then, we approximate ∇ k f n ( k ) | k α through a finitedifference approximation. As the resulting numerical gradient of the Fermi-Dirac distribution wouldbe null away from the Fermi surface, we can refine our grid in its vicinity without a substantialincrease in the computational cost, where E F is obtained by using the prescription given in Ref. 9.For those points where the numerical gradient is non-null we evaluate the components of thevelocity. Then we can replace the integral by a Riemann sum and write D ab = e π (cid:126) s (∆ k ) (cid:88) n (cid:88) α |∇ f n ( k α ) | v an ( k α ) v bn ( k α ) v n ( k α ) , (32)where ∆ k is the distance between neighbor vectors k α in our grid. III. RESULTS
We present results for the three noble metals, Ag, Au and Cu. However, we do so with moredetail for Ag, for which there is a recent accurate experimental measurement for the Drude fre-quency (cid:126) ω p = 8 . ± . Furthermore, we also consider the case of a metal under appliedstresses consisting of an isotropic and an anisotropic uniaxial strain. For the isotropic deformationwe simply modify the lattice constant a = (1 + γ ) a with respect to the equilibrium lattice constant a in the absence of stress. For the anisotropic deformation we keep the original volume of theFCC unit cell, and deform it by stretching along the z direction a z = (1 + γ ) a while shrinking italong the x and y directions, a x = a y = a (cid:107) = a / √ γ , thus converting the FCC-lattice into aBCT-lattice. We allow the expansion factor γ to take negative as well as positive values. Giventhe symmetry of these systems, we expect D ab = 0 if a (cid:54) = b , D xx = D yy ≡ D (cid:107) in the uniaxial case9nd D ab ≡ D iso δ ab in the isotropic case. Therefore, we define the Drude frequencies ω iso p = √D iso isotropic ,ω (cid:107) p = (cid:112) D (cid:107) anisotropic ,ω zp = √D z anisotropic . (33)The self-consistent ground state and the Kohn-Sham states were calculated in the DFT-LDAframework using the plane-wave ABINIT code. We used Troullier-Martins pseudopotentials that are fully separable nonlocal pseudopotentials in the Kleinman-Bylander form. We use a =4 . a = 3 . a = 4 . n picks up only the value n = 6. Finally, the number of k -points used forthe calculation was around ∼ , A. Ag
In Fig. 1 we show ω iso p ( γ ), ω (cid:107) p ( γ ) and ω zp ( γ ) as a function of the corresponding deformation γ of the unit cell. First, we point out that the three calculations agree among themselves in thelimit of no deformation γ = 0, and that the resulting value (cid:126) ω p ( γ = 0) = 9 . (cid:126) ω p = 8 . ± . . The values chosen for the compressive ( γ < γ >
0) deformations are experimentally feasible, and although small, show asizable difference for the anisotropic deformation values of ω (cid:107) p and ω zp . However, for the isotropicdeformation ω i sop does not deviate that much from its undeformed value.We notice that for very small values of γ both isotropic and anisotropic results are almost thesame, but as γ deviates away from 0, sizable changes are seen. The isotropic result ω iso p ( γ ) isrelatively flat with values between 9.2 and 9.4 eV. On the other hand, the anisotropic deformationgives an ω (cid:107) ,zp ( γ ) that decreases as a z is shortened. Also, as a z is stretched, for values of γ ≤ . ω (cid:107) p ( γ ) and ω zp ( γ ) increase and show similar values, whereas for larger values of γ ≥ . ω (cid:107) p ( γ ) continues increasing while ω zp ( γ ) reaches a maximum and starts decreasing. Our resultsshow some oscillations that are due to our approximating a surface integral through a volumeintegral represented as a sum over a discrete grid. As we deform the system, the true Fermi surfacesweeps across the grid points, giving rise to the oscillations which may be interpreted as indicativeof the accuracy of our calculation. We expect that they could be somewhat filtered away byapproximating the gradient with a higher order finite differences formula. It is worth mentioning10 . . . . . . • exp Ag shortened streched ¯ h ω p ( e V ) Deformation (%) ω iso p ω k p ω zp FIG. 1: (Color Online) Drude (cid:126) ω p vs. the deformation of the unit cell for Ag. . . . . • exp Ag shortened stretched V l + V nl only V l ω k p ¯ h ω p ( e V ) Deformation (%)
FIG. 2: (Color Online) Drude (cid:126) ω (cid:107) p vs. the deformation γ of the unit cell for Ag with and without thenonlocal potential ˆ V nl ( r , r (cid:48) ). that, as explained in Ref. 10, the experimental error in (cid:126) ω p of ∼ . H = ˆ p m + ˆ V l + ˆ V nl , (34)11 x,F compression -0.60% k y k z k x L . ×
106 m/s0 . ×
106 m/s1 . ×
106 m/s v z,F compression -0.60% k y k z k x L . ×
106 m/s0 . ×
106 m/s1 . ×
106 m/s δv F compression -0.60% k x k z k y L − . ×
106 m/s0 . ×
106 m/s
FIG. 3: (Color Online) | v x,F | (top-left), | v z,F | (top-right) and δv F = v z,F ( k x , k y , k z ) − v x,F ( k z , k y , k x )(bottom) at the Fermi Surface for Ag with an anisotropic deformation of γ = − . L point. See the text for details. where ˆ V l is the local potential characterized by the function V l ( r ), and ˆ V nl is the nonlocal potentialcharacterized by the kernel V nl ( r , r (cid:48) ). The Schr¨odinger equation reads (cid:18) − (cid:126) m ∇ + V l ( r ) (cid:19) ψ n k ( r ) + (cid:90) d r (cid:48) ˆ V nl ( r , r (cid:48) ) ψ n k ( r (cid:48) ) = E n ( k ) ψ n k ( r ) , (35)The nonlocal potential is more important for metals, where the orbitals of the conduction electronsare farther away form the nucleus, than for semiconductors, for which they are usually much closer.Therefore, we believe it is very important for our calculation to include the nonlocal contributionto the velocity v n ( k ). In our case, this was carried out using the DP code. In Fig. 2 we show ω (cid:107) p ( γ ), with and without the contribution of the nonlocal potential V nl ( r , r (cid:48) ). We see that it ismandatory to include it, otherwise the value of ω (cid:107) p would be heavily underestimated. Similar resultsare found for ω iso p and ω zp . We remark that the non-locality of the Hamiltonian arises from its beingan effective one particle effective Hamiltonian for a many-body system. As explained in Sec. II B, the numerical calculation of Eq. (32) requires the k α -points at theFS as well as the velocities evaluated at each of these points. In Fig. 3 we show the FS and theFermi velocity for an anisotropic deformation γ = − .
6% as a illustrative example. We see that12he speeds are of the order expected from the Sommerfeld model, ∼ m/s. The system issymmetric under reflections on the Cartesian planes and under an inversion around the origin. Asthe unit cell has the same lattice parameter along x and y , that v xF and v yF agree after a rotationof the system by 90 degrees in the plane, equivalent to an exchange k x ↔ k y , but they differ from v zF after a rotation by 90 degrees around the y axis due to the uniaxial strain. Fig. 3 shows that inthis case the x − z anisotropy δv F = v zF ( k x , k y , k z ) − v xF ( k z , k y , k x ) is of the order of a few percent.It is quite interesting to see that around the FS necks, centered at the L -point the speed is notuniform. B. Cu and Au
In Fig. 4 we show ω iso p ( γ ), ω (cid:107) p ( γ ) and ω zp ( γ ) as a function of the corresponding deformation γ of the unit cell for Cu and Au. As for Ag, we notice that for very small values of γ both isotropicand anisotropic results are almost the same, but as γ deviates away from 0, sizable changes areseen. The qualitative behavior of ω p for Cu and Au is very similar to that of Ag, and thus wedo not describe it again for briefness sake. The theoretical values of (cid:126) ω p ( γ = 0) = 8 .
97 eV forCu and (cid:126) ω p ( γ = 0) = 8 .
59 eV for Au for the unstrained case are close to the experimental resultsreported in Ref. 18 and shown in Table. (I). However, the experimental values, except that of Agin Ref. 10, have been reported without stating the experimental uncertainty, and there are somenoticeable differences among various experimental results, so there is definitely a need for moreprecise experimental measurements of the Drude plasma frequency for Cu and Au. However, wepoint out that our results are well within the dispersion of the available experiments.
Experimental TheoreticalThis work OthersAg 8.9 ± , 8.6 , 9.013 , 9.04 , 9.6 Cu 7.389 , 8.76 Au 7.9 , 8.55 , 8.89 , 9 , 8.951 , 9.026 (cid:126) ω p in eV for Ag, Cu, and Au in the absence of strain. C. Previous Work
We mention that in Refs. 23 and 24, the Drude plasma frequency was obtained through anexpression derived from the intraband part of the dielectric function, in the same spirit as we13 . . . . Au shortened stretched ¯ h ω p ( e V ) ω iso p ω k p ω zp . . . . . Cu shortened stretched ¯ h ω p ( e V ) Deformation (%) ω iso p ω k p ω zp FIG. 4: (Color Online) Drude (cid:126) ω p vs. the deformation of the unit cell for Au and Cu. have done here, but using a non-local longitudinal susceptibility that depends on the value of awavevector q , for which the limit q → | q | = 0 .
005 a.u. Furthermore, a fictitious finite electronic temperature T was used to smooth out the Fermi-Dirac distribution function. As we show in appendix A, takingthe actual limits q → T → (cid:126) ω p = 9 . (cid:126) ω p = 9 .
48 eV for Ag, respectively, that are not far from our values, as shown in Table(I). 14
V. CONCLUSIONS
We have shown that the well known Drude plasma frequency, ω p , should be replaced by theDrude tensor D , for which we have derived a closed theoretical expression based on quantummechanical time dependent perturbation theory. The expression for D as the integral over theFermi surface of a tensor built up from matrix elements of the velocity operator is amiable toanalytic and numerical evaluation. Using the Sommerfeld model for metals, we showed that D leads to the well known result of ω p = 4 πne /m e . We calculated D for the noble metals, and wedescribed the results for Ag in depth, since there is a recent precise result of ω p = 8 . ± . Weshowed that non-local potentials ought to be used when calculating the velocity matrix elementsin order to obtain an accurate result. In summary, the results of ω p for Cu, Ag and Au forthe undeformed unit cell with a lattice constant taken from measurements at room temperature,coincide rather well with the experimental results. We have made predictions for the noble metalssubject to an isotropic and to an uniaxial stress that could be experimentally verified. Acknowledgments
B.S.M. acknowledges the support from CONACyT through grant A1-S-9410. W.L.M. acknowl-edges the support from DGAPA-UNAM under grant IN111119. We acknowledge useful discussionswith Raksha Singla.
Appendix A: Comparison with Marini et al. Eq. (17) of Marini et al. , derived from the local limit of the longitudinal non-local suscepti-bility, reads (cid:126) ω p = lim q → πe | q | (cid:90) d k π ( f n ( k − q ) − f n ( k )) Θ ( f n ( k − q ) − f n ( k )) × (cid:12)(cid:12) (cid:104) n k | e i q · r | n k − q (cid:105) (cid:12)(cid:12) ( E n ( k − q ) − E n ( k )) , (A1)where q is the wavevector of the field. We use k · p theory to approximate this expression for small q using the following results: E n ( k + q ) = E n ( k ) + (cid:126) q · v nn ( k ) + O ( q ) | n k + q (cid:105) = e i q · r (cid:16) | n k (cid:105) + | n k (cid:105) (1) + O ( q ) (cid:17) n k (cid:105) (1) = (cid:126) (cid:88) m (cid:54) = n q · v mn ( k ) E n − E m | m k (cid:105) . (A2)Then, to first order in q , we obtain, | n k + q (cid:105) = (1 + i q · r ) (cid:16) | n k (cid:105) + | n k (cid:105) (1) + O ( q ) (cid:17) ≈ | n k (cid:105) + i q · r | n k (cid:105) + | n k (cid:105) (1) , (A3) (cid:104) n k | e i q · r | n k − q (cid:105) ≈ (cid:104) n k | (cid:16) i q · r )( | n k (cid:105) − i q · r | n k (cid:105) + | n k (cid:105) (1) (cid:17) = 1 + (cid:104) n k | n k (cid:105) (1) = 1 + (cid:126) (cid:88) m (cid:54) = n q · v mn ( k ) E n − E m (cid:104) n k | m k (cid:105) = 1 + (cid:126) (cid:88) m (cid:54) = n q · v mn ( k ) E n − E m δ nm = 1 , (A4)and E n ( k ) − E n ( k − q ) ≈ (cid:126) q · v nn ( k ) . (A5)With these results, and denoting v nn ( k ) as v n ( k ), Eq. (A1) may be written as (cid:126) ω p = lim q → πe | q | (cid:90) d k π (cid:126) q · v n ( k ) δ ( E F − E n ( k ))Θ ( f n ( k − q ) − f n ( k )) (cid:126) q · v n ( k ) ω p = e π lim q → | q | (cid:90) d k ( q · v n ( k )) δ ( E F − E n ( k )) , (A6)where the factor of 1/2 comes from the step function Θ. The expression above must be calculatedfor a given direction of q . For example, taking q along x , we get ω p = e π (cid:90) d k ( v xn ( k )) δ ( E F − E n ( k )) , (A7)equivalent to Eq. (22) for the component D xx , a the diagonal term of the Drude tensor. This isenough for isotropic systems. To get from Eq. (A6) the full Drude tensor for an anisotropic systemusing this formulation we would have to repeat the calculation for different directions of q . Also, anumerical implementation of our approach does not require a finite value of q nor a fictitious finitetemperature as used in Refs. 23 and 24. ∗ Electronic address: [email protected] P. Drude, Annalen der Physik , 566 (1900). I P. Drude, Annalen der Physik , 369 (1900). I N. W. Ashcroft and N. D. Mermin,
Solid State Physics (John Wiley & Sons, 1976). I, II A, III A G. Onida, L. Reining, and A. Rubio, Reviews of Modern Physics , 601 (2002). II S. M. Anderson, N. Tancogne-Dejean, B. S. Mendoza, and V. V´eniard, Phys. Rev. B , 075302 (2015).II C. Aversa and J. E. Sipe, , 14636 (1995-11-15), URL http://link.aps.org/doi/10.1103/PhysRevB.52.14636 . II L. H¨ormander,
The Analysis of Linear Partial Differential Operators I , vol. 256 (Springer, 1990), 2nded. II J. L. Cabellos, B. S. Mendoza, M. A. Escobar, F. Nastos, and J. E. Sipe, Phys. Rev. B , 155205 (2009).II G. A. Burdick, Phys. Rev. , 138 (1962). II B H. U. Yang, J. D’Archangel, M. L. Sundheimer, E. Tucker, G. D. Boreman, and M. B. Raschke, Phys.Rev. B , 235137 (2015). III, III A, III A, III B, III B, IV , 478 (2002); X. Gonze, et al., Zeit. Crystallogr, , 558 (2005). III N. Troullier and J. L. Martins, Phys. Rev. B , 1993 (1991). III L. Kleinman and D. M. Bylander, Phys. Rev. Lett. , 1425 (1982). III URL https://periodictable.com/Properties/A/CrystalStructure.html . III Y. Akahamaa and H. Kawamura, J. Appl. Phys. , 4767 (2004). III A V. Olevano, L. Reining, and F. Sottile, http://etsf.polytechnique.fr/Software/Ab Initio. III A O. Pulci, G. Onida, R. D. Sole, and A. J. Shkrebtii, Phys. Rev. B , 1922 (1998). III A URL . III B I. R. Hooper and J. R. Sambles, Phys. Rev. B , 165432 (2002). III B M. A. Ordal, R. J. Bell, R. W. A. Jr., L. L. Long, and M. R. Querry, Applied Optics , 4493 (1985).III B, III B, III B E. J. Zeman and G. C. Schatz, J. Phys.Chem , 634 (1987). III B, III B, III B M. G. Blaber, M. D. Arnold, and M. J. Ford, J. Phys.Chem , 3041 (2009). III B, III B A. Marini, G. Onida, and R. D. Sole, Phys. Rev. B , 115101 (2002). III B, III C, III C, A A. Marini, G. Onida, and R. D. Sole, Phys. Rev. B , 195125 (2001). III B, III C, III C, A, A M. Kreiter, S. Mittler, W. Knoll, and J. Sambles, Phys. Rev. B , 125415 (2002). III B S. Berciaud, L. Cognet, P. Tamarat, and B. Lounis, Nano Letters , 515 (2005). III B N. K. Grady, N. J. Halas, and P. Nordlander, Chem. Phys. Lett. , 167 (2004). III B, 167 (2004). III B