Theory of non-retarded ballistic surface plasma waves in metal films

TL;DR

This paper develops a comprehensive theory for surface plasma waves in metal films, revealing a universal self-amplification mechanism unaffected by electronic collision rates, with implications for experimental phenomena like Ferrel radiation.

Contribution

It introduces a novel, universal self-amplification channel for surface plasma waves in metal films, independent of collision rates, and analyzes effects of film thickness and surface roughness.

Findings

Discovery of a universal self-amplification channel for plasma waves

Independence of amplification properties from electronic collision rate

Implications for Ferrel radiation and experimental observations

Abstract

We present a theory of surface plasma waves in metal films with arbitrary electronic collision rate $\tau$. Both \textit{tangential} and \textit{normal} modes are investigated. A universal self-amplification channel for these waves is established as a result of the unique interplay between ballistic electronic motions and boundary effects. The channel is shown to be protected by a general principle and its properties independent of $\tau$. The effects of film thickness and surface roughness are also calculated. Experimental implications, such as Ferrel radiation, are discussed.

Equations

$$\left({\partial }_t+\frac{1}{\tau }\right)\rho (\mathbf{x},t)+{\partial }_{\mathbf{x}}\cdot \mathbf{J}(\mathbf{x},t)=S(\mathbf{x},t),$$

$$S(\mathbf{x},t)=J_z({\mathbf{x}}_d,t)\delta (z-d)-J_z({\mathbf{x}}_0,t)\delta (z),$$

$$-i\bar{\omega }\rho (z)+\nabla \cdot \mathbf{J}(z)=S(z),$$

$$S(z)=J_z(d)\delta (z-d)-J_z(0)\delta (z).$$

$$\varphi (z)=\frac{2\pi }{k}\int d{z}^{\prime }{e}^{-k|z-z^{\prime }|}\rho (z^{\prime }).$$

$$\rho (z)=\sum _{n=0}^{\infty }{\rho }_n\cos (q_{n}z),q_n=\frac{\pi n}{d}.$$

$${\rho }_n=\frac{1}{d_n}{\int }_0^{d}dz\rho (z)\cos (q_{n}z),d_n=\frac{d}{2-{\delta }_{n,0}},$$

$$\rho (z)={\rho }_{+}(z)+{\rho }_{-}(z),$$

$$\varphi (z)=\sum _{n=0}^{\infty }\frac{2\pi {\rho }_n}{k^2+q_n^2}\left[2\cos (q_{n}z)-e^{-kz}-(-1{)}^n{e}^{-k(d-z)}\right],$$

$$\frac{\partial g(\mathbf{v},z)}{\partial z}+{\lambda }^{-1}g(\mathbf{v},z)+e{f}_0^{\prime }(\epsilon )\frac{\mathbf{v}\cdot \mathbf{E}(z)}{v_z}=0,$$

$$g(\mathbf{v},z)=g_{bulk}(\mathbf{v},z)+g_{surface}(\mathbf{v},z),$$

$$g_{bulk}(\mathbf{v},z)=-e{f}_0^{\prime }\sum _{n=-\infty }^{\infty }\frac{2\pi {\rho }_n}{k^2+q_n^2}\frac{k{v}_x+q_n{v}_z}{\bar{\omega }-(k{v}_x+q_n{v}_z)}{e}^{i{q}_{n}z},$$

$$g_{D,surface}(\mathbf{v},z)=g_{D,surface}^{(1)}(\mathbf{v},z)+g_{D,surface}^{(2)}(\mathbf{v},z),$$

$$g_{D,surface}^{(1)}(\mathbf{v},z)=-e{f}_0^{\prime }\sum _{n=0}^{\infty }\frac{2\pi {\rho }_n}{k^2+q_n^2}\frac{k(v_z-i{v}_x)}{k{v}_z+i\tilde{\omega }}{e}^{-kz},$$

$$g_{D,surface}^{(2)}(\mathbf{v},z)=-e{f}_0^{\prime }\sum _{n=0}^{\infty }\frac{2\pi {\rho }_n(-1{)}^n}{k^2+q_n^2}\frac{k(v_z+i{v}_x)}{k{v}_z-i\tilde{\omega }}{e}^{-k(d-z)}.$$

$$g_D(\mathbf{v},z)=g_{bulk}(\mathbf{v},z)+g_{D,surface}(\mathbf{v},z),$$

$$g_B(\mathbf{v},z)\text{mathchar}58=g_{B,surface}(\mathbf{v},z).$$

$$g(\mathbf{v},z)=g_D(\mathbf{v},z)+g_B(\mathbf{v},z)$$

$$g_{B,surface}(\mathbf{v},z)=g_{B,surface}^{(1)}(\mathbf{v},z)+g_{B,surface}^{(2)}(\mathbf{v},z),$$

$$g_{B,surface}^{(1)}(\mathbf{v},z)=\Theta (v_z)e^{i\frac{\tilde{\omega }z}{v_z}}\left[g_{B,emg}^{(1)}(\mathbf{v})+p_1{g}_{B,ref}^{(1)}(\mathbf{v})\right],g_{B,surface}^{(2)}(\mathbf{v},z)=\Theta (-v_z)e^{i\frac{\tilde{\omega }(z-d)}{v_z}}\left[g_{B,emg}^{(2)}(\mathbf{v})+p_2{g}_{B,ref}^{(2)}(\mathbf{v})\right],$$

$$g_{B,emg}^{(1/2)}(\mathbf{v})=e{f}_0^{\prime }\sum _{n=0}^{\infty }\frac{2\pi {\rho }_n⟨1|(-1{)}^n⟩}{k^2+q_n^2}\left[\frac{2(\tilde{\omega }k{v}_x)+q_n^2{v}_z^2}{{\tilde{\omega }}^2-q_n^2{v}_z^2}+\frac{k(v_z\mp i{v}_x)}{k{v}_z\pm i\tilde{\omega }}+(-1{)}^n{e}^{-kd}\frac{k(v_z\pm i{v}_x)}{k{v}_z\mp i\tilde{\omega }}\right],$$

$$g_{B,ref}^{(1)}(\mathbf{v})=e{f}_0^{\prime }\sum _{n=0}^{\infty }\frac{2\pi {\rho }_n}{k^2+q_n^2}\left[\left(e^{i\frac{\tilde{\omega }d}{v_z}}(-1{)}^n-1\right)\frac{2(\tilde{\omega }k{v}_x-q_n^2{v}_z^2)}{{\tilde{\omega }}^2-q_n^2{v}_z^2}+\left(1-e^{\left(\frac{i\tilde{\omega }}{v_z}-k\right)d}\right)\frac{k(v_z-i{v}_x)}{k{v}_z-i\tilde{\omega }}+(-1{)}^n\left(e^{-kd}-e^{i\frac{\tilde{\omega }d}{v_z}}\right)\frac{k(v_z+i{v}_x)}{k{v}_z+i\tilde{\omega }}\right],$$

$$g_{B,ref}^{(2)}(\mathbf{v})=e{f}_0^{\prime }\sum _{n=0}^{\infty }\frac{2\pi {\rho }_n(-1{)}^n}{k^2+q_n^2}\left[\left(e^{-i\frac{\tilde{\omega }d}{v_z}}(-1{)}^n-1\right)\frac{2(\tilde{\omega }k{v}_x-q_n^2{v}_z^2)}{{\tilde{\omega }}^2-q_n^2{v}_z^2}+\left(1-e^{-\left(\frac{i\tilde{\omega }}{v_z}+k\right)d}\right)\frac{k(v_z+i{v}_x)}{k{v}_z+i\tilde{\omega }}+(-1{)}^n\left(e^{-kd}-e^{-i\frac{\tilde{\omega }d}{v_z}}\right)\frac{k(v_z-i{v}_x)}{k{v}_z-i\tilde{\omega }}\right],$$

$${\mathbf{J}}_{D/B}(z)=(m/2\pi \hslash {)}^3\int d^3\mathbf{v}e\mathbf{v}{g}_{D/B}(\mathbf{v},z)$$

$${\mathbf{J}}_{D,surface}=\frac{i}{\bar{\omega }}\frac{{\omega }_p^2}{4\pi }{\mathbf{E}}^{surface}(z),$$

$${\mathbf{J}}_{bulk}(z)=\frac{i}{\bar{\omega }}\frac{{\omega }_p^2}{4\pi }{\mathbf{E}}^{bulk}(z)+{\mathbf{J}}^{\prime }(z),$$

$${\mathbf{J}}^{\prime }(z)=\sum _{n=-\infty }^{\infty }\frac{2\pi {\rho }_n{e}^{i{q}_{n}z}}{k^2+q_n^2}\mathbf{F}(k,q_n;\bar{\omega }).$$

$$\mathbf{F}(k,q;\bar{\omega })={\left(\frac{m}{2\pi \hslash }\right)}^3\int d^3\mathbf{v}(-e^2{f}_0^{\prime })\mathbf{v}\sum _{l=2}^{\infty }{\left(\frac{k{v}_x+q{v}_z}{\bar{\omega }}\right)}^l.$$

$${\mathbf{J}}_D(z)=\frac{i}{\bar{\omega }}\frac{{\omega }_p^2}{4\pi }\mathbf{E}(z)+{\mathbf{J}}^{\prime }(z).$$

$$\nabla \cdot {\mathbf{J}}_D(z)=\frac{i}{\bar{\omega }}\sum _{n=0}^{\infty }{\Omega }^2(k,q_n;\bar{\omega }){\rho }_n\cos (q_{n}z),$$

Full text

I Introduction

Surface plasma waves (SPWs) ritchie1957 ; ferrell1958 ; raether1988 are fascinating to a wide spectrum of scientists not only for their fundamental physical properties pitarke2007 ; feibelman1982 ; pendry1975 but also their promising potential int in a myriad of applications, including microscopy benno1988 , sensing hoa2007 and nano-optics zayats2005 ; maier2007 ; ozbay2006 ; mark2007 ; barnes2003 ; ebbesen2008 as well as information processing tame2013 . Being charge density waves highly localized about the interface between a metal and a dielectric, SPWs strongly interact and form a bound entity with light that might render an atomic resolution of molecular dynamics nie1997 . Nowadays SPWs are pivotal in nano-optics.

The standard theory of SPWs was delivered shortly after the pioneering work ritchie1957 by Ritchie in 1957 and has since been comprehensively discoursed in many textbooks and review articles raether1988 ; pitarke2007 ; maier2007 ; sarid2010 . In this theory, the electrical properties of a metal are prescribed with a dielectric function $\epsilon$. To analytically treat $\epsilon$, the simple Drude model or the slightly more involved hydrodynamic model is often invoked pitarke2007 ; harris1971 ; fetter1986 ; pendry2013 ; pendry2014 ; pendry2015 ; schnitzer2016 . For either model to be valid, electronic collisions in the metal must be sufficiently frequent so that the electronic mean free path, $l_{0}=v_{F}\tau$, where $v_{F}$ is the Fermi velocity and $\tau$ the thermal charge relaxation time, is much shorter than the SPW wavelength or the typical length of the system ziman1960 ; pines ; abrikosov . The general case with arbitrary $\tau$, especially the collision-less limit, where $\tau\rightarrow\infty$, defies these models and has yet to be entertained. Other models based on ab initio quantum mechanical computations pitarke2007 are helpful in understanding the complexity of real materials but falls short in providing an intuitive and systematic picture of SPWs underpinned by electrons experiencing less frequent collisions.

The purpose of this paper is to furnish a comprehensive theory for SPWs of ballistically moving electrons. Ballistic SPWs are not only interesting in themselves but could have ramified applications in plasmonics and other arenas. Recently deng2017a ; deng2017b , we considered ballistic SPWs in semi-infinite metals. We showed that such waves are intrinsically unstable and possess a universal self-amplification channel that exists irrespective of the value of $\tau$. This result was initially established by examining the charge dynamics deng2017a in the system and later corroborated by an energy conversion analysis deng2017b in the waves. In the present work, we study ballistic SPWs in metal films, which possess two surfaces and are experimentally more realistic and interesting.

In the next section, we specify the system under consideration and state our main results. Some preliminary remarks are made on their experimental implications. In Sec. III, the theory in support of the results is systematically presented, followed by a complementary energy conversion analysis in Sec. IV. We discuss the results and conclude the paper in Sec. V. Some calculations of technical interest are displayed in the appendices A, B and C.

II Results

System. We consider ballistic SPWs in a metal film surrounded by vacuum. By the so-called jellium model pines ; abrikosov , the metal is described as a free electron gas embedded in a static background of homogeneously distributed positive charges. This description is valid if the length scale in question is much longer than the microscopic lattice constant and inter-band transitions are negligible. The kinetic energy of electrons is $\varepsilon(\mathbf{v})=\frac{1}{2}m\mathbf{v}^{2}$, where $m$ and $\mathbf{v}$ denote the mass and velocity of the electrons, respectively. The film resides in the region $0\leq z\leq d$ with two surfaces located at $z=0$ and $z=d$, respectively. The surfaces are treated as geometric planes of a hard wall type and they strictly prevent electrons from leaking out of the metal. To simplify our analysis, the surfaces are assumed with identical properties so that the system is symmetric about the mid-plane $z=d/2$. To avoid quantum size effects, we assume $d\gg\hbar/mv_{F}$, where $\hbar$ is the reduced Planck constant. Throughout we write $\mathbf{x}=(\mathbf{r},z)$ and reserve $\mathbf{r}=(x,y)$ for planar components while let $t$ be the time. We neglect retardation effects in total fetter1986 ; deng2015 .

Results. With two surfaces, a film possesses two branches of SPWs, which at large $d$ degrade into those for two semi-infinite metals. Reflection symmetry about the mid-plane requires the corresponding charge densities to bear a definite sign under the reflection. The branch whose charge density is invariant under the reflection is called symmetric while the one whose charge density changes sign under reflection is called anti-symmetric. In the literature, the symmetric and anti-symmetric SPWs are also designated as tangential and normal oscillations, respectively. Profiles of the charge densities for symmetric and anti-symmetric SPWs are mapped in Fig. 1 (a) and (b), respectively, together with the electric field accompanying them.

We find that the SPW frequency $\omega^{\pm}_{s}$ is significantly (as much as $30$%) higher than $\omega^{\pm}_{s0}=\left(\omega_{p}/\sqrt{2}\right)\sqrt{1\pm e^{-kd}}$ which would be obtained by the hydrodynamic/Drude theory. Here the plus (minus) sign is affixed and refers to symmetric (anti-symmetric) modes, $\omega_{p}$ denotes the characteristic plasma frequency of the metal and $k$ is the SPW wavenumber. The dependences of $\omega^{\pm}_{s}$ on $k$, $d$ and surface scattering – the effects of which could be summarized in the Fuchs parameter $p$ in the simplest possible scattering picture, are displayed in the upper panels of Fig. 2 (a), (b) and (c), respectively. The great contrast between $\omega^{\pm}_{s}$ and $\omega^{\pm}_{s0}$ would be ideal for experimentally verifying our theory. Unfortunately, in the most commonly experimented materials, such as noble metals, due to pronounced inter-band transitions there is no simple relation between $\omega_{p}$ and $\omega^{\pm}_{s0}$.

More interestingly, we reveal a universal self-amplification channel for SPWs irrespective of their symmetry. Namely, we find that the net amplification rate of SPWs can be generally written as $\gamma^{\pm}=\gamma^{\pm}_{0}-\tau^{-1}$, where $\gamma_{0}$ is warranted to be non-negative by a general principle and independent of $\tau$. In the conventional theory, $\gamma^{\pm}_{0}$ vanishes identically and amplification would be impossible without extrinsic energy supply bergman2003 ; seidel2005 ; leon2008 ; leon2010 ; yu2011 ; pierre2012 ; fedyanin2012 ; cohen2013 . The dependences of $\gamma^{\pm}_{0}$ on $k$, $d$ and $p$ are shown in the lower panels of Fig. 2 (a), (b) and (c), respectively, where we observe that (1) $\gamma^{\pm}_{0}$ is generally a sizable fraction (as much as $\sim 10$%) of $\omega_{p}$, (2) it increases as $k$ increases, i.e. higher amplification obtains for shorter wavelengths and (3) it increases as $p$ increases, i.e. smooth surfaces produce higher amplification than rough surfaces. We also see that $\gamma^{+}_{0}$ is more sensitive to film thickness than $\gamma^{-}_{0}$.

Additionally, we show that the electrical current density in the system can be split into two disparate components, which we call $\mathbf{J}_{D}$ and $\mathbf{J}_{B}$, respectively. An example of their profiles is exhibited in Fig. 3 (a) and (b) respectively for the symmetric and anti-symmetric modes. What critically sets them apart rests with their distinct relations with the electric field $\mathbf{E}$ present in the system. $\mathbf{J}_{D}$ responds to $\mathbf{E}$ as if the system had no surfaces and is therefore primarily a bulk property. As such, it can also be satisfactorily captured by the hydrodynamic/Drude model. For this reason, we designate it a diffusive component, regardless of the value of $\tau$. On the contrary, $\mathbf{J}_{B}$ represents genuine surface effects and would totally disappear were the surfaces absent. In particular, it synthesizes the effects ensuing from the fact that the system is not translationally invariant along the direction normal to the surfaces. These effects are completely beyond the hydrodynamic/Drude model but well within the scope of Boltzmann’s approach, which is employed in our theory to be expounded in the next section. We thus designate $\mathbf{J}_{B}$ as a surface-ballistic component.

Finally, we find that the self-amplification channel is a direct consequence of $\mathbf{J}_{B}$. Indeed, were not for $\mathbf{J}_{B}$, SPWs would behave in accord with the hydrodynamic/Drude model. This is already clear from the orientations of $\mathbf{J}_{D/B}$ relative to $\mathbf{E}$. As seen in Fig. 3, $\mathbf{J}_{D}$ points at right angles with $\mathbf{E}$ almost locally, whereas $\mathbf{J}_{B}$ flows normal to the surface paying little regard to $\mathbf{E}$. Therefore, $\mathbf{E}$ does no work on $\mathbf{J}_{D}$ on average while, as shown in Sec. IV, it does a negative amount of work on $\mathbf{J}_{B}$, thereby imparting energy from the electrons to SPWs and destabilizing the Fermi sea.

Remarks. Experimentally verifying the self-amplification channel and the theory in general would be of considerable interest, as it would drastically change the way we conceive and utilize SPWs and renew our interest in surface science in a broad sense. The self-amplification channel could manifest itself for instance in the temperature dependence of various spectra, e.g. electron loss spectra. We discuss this aspect in Sec. V. Here we mainly concern ourselves with the experimental implications of the surface-ballistic current $\mathbf{J}_{B}$.

Being an integral part of the electrical responses of metals, $\mathbf{J}_{B}$ is expected to play a role in virtually every phenomena where surface is not negligible. Examples include electron energy losses, reflectance and van der Waals forces. Unlike $\mathbf{J}_{D}$, which does not reflect surface scattering effects, $\mathbf{J}_{B}$ is surface-specific via the Fuchs parameter. Moreover, they differ in phase by $\sim\pi/2$. To be specific, let us consider the Ferrel radiation ferrell1958 . Ferrel predicted that anti-symmetric SPWs in thin films would radiate in a characteristic pattern. Some experiments even claimed to have observed this radiation brown1960 ; arakawa1964 ; donohue1986 . Ferrel considered only $\mathbf{J}_{D}$. Following him, we find that including $\mathbf{J}_{B}$ could boost the radiation power by a factor $\sim 1+(3/2\pi)^{2}(1+2p)^{2}$. Though a crude estimate, it does imply that surface properties could be utilized to tune the radiation. In this paper, we focus on the fundamental theory of ballistic SPWs. A systematic treatment of Ferrel radiation will be published elsewhere.

As aforementioned, a major obstacle in experimentally studying the theory lies with inter-band transitions, which have been neglected in our theory. A detailed discussion of their effects is presented in Sec. V.

III Theory

This section is devoted to a thorough exposition of the theory. We begin with a discussion of the equation of continuity in the presence of surfaces. Thence we proceed to Boltzmann’s approach and analyze how to handle surface effects in this approach. The electronic distribution functions, obtained by solving Boltzmann’s equation, are discussed in detail. The electrical current densities are then calculated and the exact equation of motion for the charge density is established. Solutions to the equation are discussed and the properties of SPWs are analyzed. Various limits are presented and connections are made with the hydrodynamic/Drude models.

III.1 Equation of Continuity

The starting point of our theory is the equation of continuity, $\left(\partial_{t}+1/\tau\right)\rho(\mathbf{x},t)+\partial_{\mathbf{x}}\cdot\mathbf{j}(\mathbf{x},t)=0$, which relates the charge density $\rho(\mathbf{x},t)$ and the current density $\mathbf{j}(\mathbf{x},t)$ in a universal manner. Here $\mathbf{j}(\mathbf{x},t)$ arises in the presence of an electric field $\mathbf{E}(\mathbf{x},t)$ and the damping term $-\rho(\mathbf{x},t)/\tau$ is included to account for the thermal currents due to electronic collisions that would drive the system toward thermodynamic equilibrium. In the jellium model, $\rho(\mathbf{x},t)$ appears when the electron density is perturbed away from its equilibrium value $n_{0}$.

As the surfaces strictly prevent electrons from escaping the metal, we may write $\mathbf{j}(\mathbf{x},t)=\left[\Theta(z)-\Theta(z-d)\right]\mathbf{J}(\mathbf{x},t)$, where $\Theta(z)$ is the Heaviside step function. In doing this, we have embodied the surfaces as hard walls and considered the fact that $\mathbf{J}(\mathbf{x},t)$ may not vanish even in the immediate neighborhood of the surfaces – as is obviously the case with Drude model. With this prescription, the equation of continuity can be rewritten

[TABLE]

where the effective source term

[TABLE]

results directly from the presence of the surfaces. Here $\mathbf{x}_{0}=(\mathbf{r},0)$ and $\mathbf{x}_{d}=(\mathbf{r},d)$ denote points on the surface at $z=0$ and those on that at $z=d$, respectively. Physically, $S(\mathbf{x},t)$ corresponds to the scenario that charges must pile up on the surfaces if they do not come to a halt before they reach them.

Without loss of generality we seek fields in this form: $\rho(\mathbf{x},t)=$ Re $\left[\rho(z)e^{i(kx-\omega t)}\right]$ and $\mathbf{J}(\mathbf{x},t)=$ Re $\left[\mathbf{J}(z)e^{i(kx-\omega t)}\right]$. Similarly, for the electric field $\mathbf{E}(\mathbf{x},t)=$ Re $\left[\mathbf{E}(z)e^{i(kx-\omega t)}\right]$ and the electrostatic potential $\phi(\mathbf{x},t)=$ Re $\left[\phi(z)e^{i(kx-\omega t)}\right]$. In these expressions, Re/Im takes the real/imaginary part of a quantity, $k\geq 0$ is a wavenumber and $\omega$ is the eigen-frequency to be determined. Equation (1) becomes

[TABLE]

where $\bar{\omega}=\omega+i/\tau$, $\mathbf{\nabla}=(ik,\partial_{y},\partial_{z})$ and

[TABLE]

Equation (3) will serve as the equation of motion for $\rho(z)$ when supplemented with additional relations to be formulated between $\mathbf{J}(z)$ and $\rho(z)$ in what follows.

III.2 The Law of Electrostatics

If the SPW phase velocity is much smaller than the speed of light $c$ in vacuum, i.e. $k>k_{0}$, where $k_{0}$ is the wavenumber of light at the SPW frequency, the system will be in the non-retarded regime deng2015 and we can relate $\phi(\mathbf{x},t)$ and $\rho(\mathbf{x},t)$ by the laws of electrostatics. Without external charges, we have deng2015

[TABLE]

Instead of $\rho(z)$, we directly work with its Fourier components. Generically, we may write

[TABLE]

The components are given by

[TABLE]

where $\delta_{m,n}$ denotes the Kroneker symbol.

As the surfaces of the film are assumed identical, the system is invariant under reflection about its mid-plane. This symmetry makes it useful to write $\rho(z)$ as a superposition of a symmetric mode $\rho_{+}(z)$ and an anti-symmetric mode $\rho_{-}(z)$. Namely,

[TABLE]

where $\rho_{+}(z)$ includes all the terms with even $n$ whereas $\rho_{-}(z)$ those with odd $n$. As such, $\rho_{+}(0)=\rho_{+}(d)$ and $\rho_{-}(0)=-~{}\rho_{-}(d)$. Due to the symmetry $\rho_{+}(z)$ and $\rho_{-}(z)$ will be shown to be strictly decoupled. We impose on $q_{n}$ a cutoff $q_{c}$ of the order of a reciprocal lattice constant; otherwise, the jellium model would cease to be valid. Obviously, $q_{c}\sim k_{F}$, where $k_{F}$ is the Fermi wavenumber of the electrons in the metal.

In terms of $\rho_{n}$, we can rewrite

[TABLE]

The electric field, $\mathbf{E}(z)=-~{}\mathbf{\nabla}\phi(z)$, can then be obtained straightforwardly. In equation (6) the exponentials, $e^{-kz}$ and $e^{-k(d-z)}$, would all vanish if the surfaces were sent to infinity. We may then write $\phi(z)=\phi^{bulk}(z)+\phi^{surface}(z)$, where $\phi^{surface}(z)$ includes the contributions from all the exponentials while $\phi^{bulk}(z)$ contains the remaining contributions. Accordingly, $\mathbf{E}(z)=\mathbf{E}^{bulk}(z)+\mathbf{E}^{surface}(z)$. Such a partition proves useful in analyzing surface specific effects.

III.3 Electronic Distribution Function

The electric field $\mathbf{E}(\mathbf{x},t)$ drives an electrical current $\mathbf{J}(\mathbf{x},t)$. We employ Boltzmann’s equation, which is valid as long as inter-band transitions are negligible, to calculate this current. Including the transitions in our formalism is straightforward but will be skipped here. Surfaces scatter electrons. On the microscopic level, one can in principle introduce a surface potential $\phi_{s}(\mathbf{x})$ in Boltzmann’s equation to produce such scattering. The corresponding surface field $\mathbf{E}_{s}(\mathbf{x})=-\partial_{\mathbf{x}}\phi_{s}(\mathbf{x})$ should be peaked on the surfaces and may have an infinitesimal spread complying with the hard-wall picture of surfaces. However, as $\phi_{s}(\mathbf{x})$ can hardly be known and varies from one sample to another, this method is impractical and futile.

Alternatively surface scattering effects can be dealt with using boundary conditions. This is possible because $\mathbf{E}_{s}(\mathbf{x})$ acts only within the immediate neighborhoods of the surfaces; In the bulk of the sample, the electronic distribution function $f(\mathbf{x},\mathbf{v},t)$ sought as solutions to Boltzmann’s equation can be specified up to some parameters, which summarize the effects of – while without actually knowing – $\phi_{s}(\mathbf{x})$. With translational symmetry along the surfaces, only one such parameter, i.e. the so-called Fuchs parameter $p$, is needed in the simplest model. Physically, $p$ measures the probability that an electron is bounced back when impinging upon the surface. We write $f(\mathbf{x},\mathbf{v},t)=f_{0}(\varepsilon(\mathbf{v}))+g(\mathbf{x},\mathbf{v},t)$, where $f_{0}(\varepsilon)$ denotes the Fermi-Dirac distribution and $g(\mathbf{x},\mathbf{v},t)$ represents the non-equilibrium part due to the presence of $\mathbf{E}(\mathbf{x},t)$. The current density can then be calculated by $\mathbf{J}(\mathbf{x},t)=(m/2\pi\hbar)^{2}\int d^{3}\mathbf{v}~{}e\mathbf{v}~{}g(\mathbf{x},\mathbf{v},t)$, where $e$ denotes the charge of an electron. It is worth pointing out that, as $g(\mathbf{x},\mathbf{v},t)$ is a distribution for the bulk, the actual charge density is not given by $\tilde{\rho}(\mathbf{x},t)=(m/2\pi\hbar)^{2}\int d^{3}\mathbf{v}~{}e~{}g(\mathbf{x},\mathbf{v},t)$, i.e. $\rho(\mathbf{x},t)\neq\tilde{\rho}(\mathbf{x},t)$. Actually, $\tilde{\rho}(\mathbf{x},t)$ satisfies $(\partial_{t}+1/\tau)\tilde{\rho}(\mathbf{x},t)+\partial_{\mathbf{x}}\cdot\mathbf{J}(\mathbf{x},t)=0$ rather than Eq. (1). By comparison, one sees that what is missing from $\tilde{\rho}(\mathbf{x},t)$ is the charges localized on the surface.

As before we write $g(\mathbf{x},\mathbf{v},t)=$ Re $\left[g(\mathbf{v},z)e^{i(kx-\omega t)}\right]$. For linear responses, Boltzmann’s equation can be written

[TABLE]

where $\lambda=iv_{z}/\tilde{\omega}$ with $\tilde{\omega}=\bar{\omega}-kv_{x}$ and $f^{\prime}_{0}(\varepsilon)=\partial f_{0}/\partial\varepsilon(\mathbf{v})$. In this equation, the velocity $\mathbf{v}$ is more of a parameter than an argument and can be used to tag electron beams. It is straightforward to solve the equation under appropriate boundary conditions seeAppendixB. We divide $g(\mathbf{v},z)$ into a bulk and a surface term, i.e.

[TABLE]

where the bulk term would exist even in the absence of surfaces whereas the surface term would not. Using Eq. (6) for $\mathbf{E}(z)$, we obtain

[TABLE]

where we have defined $\rho_{n<0}\mathrel{\mathop{\mathchar 58\relax}}=\rho_{-n}$. For large $d$ equation (8) converges to the distribution function of a boundless system for either the symmetric mode or the anti-symmetric mode. It is notable that $g_{bulk}(\mathbf{v},z)$ bears a single form for all electrons regardless of their velocities.

As for $g_{surface}(\mathbf{v},z)$, we find it with a subtle structure: it can be written as a sum of two contributions, one of which, $g_{D,surface}(\mathbf{v},z)$, has a single form for all electrons irrespective of their velocities while the other, $g_{B,surface}(\mathbf{v},z)$, does not. Explicitly, we find

[TABLE]

where

[TABLE]

and

[TABLE]

originate from the surfaces at $z=0$ and $z=d$, respectively.

We may combine $g_{bulk}(\mathbf{v})$ and $g_{D,surface}(\mathbf{v},z)$ in a single term,

[TABLE]

in order to separate them from

[TABLE]

The subscripts, $D$ and $B$, refer to ’diffusive’ and ’surface-ballistic’, respectively. In so doing, we have decomposed

[TABLE]

in a diffusive and a surface-ballistic component. It is underlined that $g_{B}(\mathbf{v},z)$ arises only when the surfaces are present. For boundless systems without surfaces, it does not exist even if the electronic motions are totally ballistic, i.e. $\tau\rightarrow\infty$. In other words, $g_{B}(\mathbf{v},z)$ represents genuine surface effects. It may be interpreted as a contribution from electrons which experience the electric field only on the surfaces and propagate freely in the body. Its expressions are given in what follows.

Electrons in the film can bounce back and forth between its surfaces. Each bounce gives a factor $p_{1}p_{2}e^{2i\tilde{\omega}d/\mathinner{\!\left\lvert v_{z}\right\rvert}}$, whose magnitude is generally smaller than unity (see Appendix B). Here $p_{1}$ and $p_{2}$ are the Fuchs parameters for the surfaces at $z=0$ and $z=d$, respectively. Consequently, we neglect multiple bounces, which allows to write

[TABLE]

where $g^{(1)}_{B,surface}(\mathbf{v},z)$ and $g^{(2)}_{B,surface}(\mathbf{v},z)$ originate from the surfaces at $z=0$ and $z=d$, respectively. They are given by

[TABLE]

where $g^{(1/2)}_{B,emg}(\mathbf{v},z)$ is contributed by electrons that directly emerge from the surface at $z=0/d$ while $g^{(1/2)}_{B,ref}(\mathbf{v},z)$ by reflected electrons and hence proportional to $p_{1/2}$. In what follows we take $p_{1}=p_{2}=p$. The expressions of $g^{(1/2)}_{B,emg/ref}(\mathbf{v})$ are involved but with a recognizable structure:

[TABLE]

where the symbol $\langle 1|(-1)^{n}\rangle$ returns $1$ and $(-1)^{n}$ for $g^{(1)}_{B,emg}$ and $g^{(2)}_{B,emg}$, respectively. In addition, we have

[TABLE]

and

[TABLE]

Positiveness of $\mbox{Im}(\bar{\omega})$. What sets $g_{B}(\mathbf{v},z)$ apart from its diffusive counterpart rests with its disparate $z$ dependence. Let us take the contribution originating from the surface at $z=0$ for example. Here $g^{(1)}_{B,surface}(\mathbf{v},z)\propto e^{i\tilde{\omega}z/v_{z}}\propto e^{-\mbox{Im}(\bar{\omega})z/v_{z}}$, where $v_{z}\geq 0$. Unless Im$(\bar{\omega})\geq 0$, this expression would diverge for small $v_{z}$. As such, we may conclude that Im$(\bar{\omega})\geq 0$, a result to be confirmed in what follows by specific calculations. In Appendix B, we frame this result as a consequence of the causality principle: out-going electrons are determined by in-coming ones; not otherwise.

III.4 Current Densities

We are now prepared to discuss the behaviors of the current density, which is written $\mathbf{J}(z)=\mathbf{J}_{D}(z)+\mathbf{J}_{B}(z)$, where

[TABLE]

is the diffusive/surface-ballistic component of $\mathbf{J}(z)$. The equation of motion for $\rho(z)$ follows upon inserting $\mathbf{J}(z)$ in Eq. (3). In our calculations, the zero temperature is assumed whenever a concrete form of $f_{0}(\varepsilon)$ is required, though generalization to finite temperatures is straightforward.

III.4.1 Diffusive current density

Since $g_{D}(\mathbf{v},z)$ consists of a bulk and a surface component, we accordingly write $\mathbf{J}_{D}(z)=\mathbf{J}_{bulk}(z)+\mathbf{J}_{D,surface}(z)$, where $\mathbf{J}_{bulk}(z)$ and $\mathbf{J}_{D,surface}(z)$ arise from $g_{bulk}(\mathbf{v},z)$ and $g_{D,surface}(\mathbf{v},z)$, respectively. By straightforward manipulation, one may show that $\mathbf{J}_{D,surface}(z)\propto\mathbf{E}^{surface}(z)$. To the lowest order in $kv_{F}/\omega_{p}$, where $\omega_{p}=\sqrt{4\pi n_{0}e^{2}/m}$ is the characteristic plasma frequency of the metal, we have

[TABLE]

where the pre-factor heading $\mathbf{E}^{surface}(z)$ is recognized as the Drude conductivity. In addition, we find

[TABLE]

where

[TABLE]

signifies non-local electrical responses that would engender dispersive plasma waves. In the expression

[TABLE]

Only terms with odd $l$ contribute in the series. Note that the normal component of $\mathbf{J}^{\prime}(z)$ vanishes identically at all surfaces, i.e. $J^{\prime}_{z}(0)=J^{\prime}_{z}(d)=0.$

Piecing everything together we obtain

[TABLE]

As in the hydrodynamic/Drude model, which is valid only for diffusive electronic motions, the relation between $\mathbf{J}_{D}(z)$ and $\mathbf{E}(z)$ assumes the form of a generalized Ohm’s law. This is why we consider $\mathbf{J}_{D}(z)$ a diffusive component, irrespective of the value of $\tau$. Its divergence is easily found to be

[TABLE]

where, with $\mathbf{k}\mathrel{\mathop{\mathchar 58\relax}}=(k,q)$,

[TABLE]

Fourier transforming Eq. (18) yields

[TABLE]

We will show that $\Omega(k,q;\bar{\omega})$ is intimately related to the properties of bulk plasma waves. As expected, $\Omega(k,q;\bar{\omega})$ only depends on the length of $\mathbf{k}$, not its direction. This becomes evident by writing $kv_{x}+qv_{z}=\mathbf{k}\cdot\mathbf{v}$ in Eq. (17). The first non-vanishing contribution to $\Omega(k,q;\bar{\omega})$ comes from the term $l=1$ in the series in $\mathbf{F}(k,q;\bar{\omega})$. Retaining only this term, we get

[TABLE]

Upon replacing $\bar{\omega}$ with $\omega_{p}$, one immediately revisits the dispersion relation for bulk waves, which could also be reached through the hydrodynamic model. In the Drude model, the dispersion is totally neglected.

It is noted that $\Omega(k,q;\bar{\omega})$ generally possesses an imaginary part. In case Im$(\bar{\omega})$ is vanishingly small, the imaginary part arises from a pole, located at $\bar{\omega}=kv_{x}+qv_{z}$, in the integrand in Eq. (17), giving rise to Landau damping in bulk waves and SPWs. In our numerical computation of $\bar{\omega}$, Landau damping will be automatically included.

III.4.2 Surface-ballistic current density

Separating the contributions of emerging electrons from that of reflected electrons, we write $\mathbf{J}_{B}(z)=\mathbf{J}_{B,emg}(z)+p\mathbf{J}_{B,ref}(z)$. Explicitly, we find

[TABLE]

where we have defined $\mathbf{J}_{B,emg/ref}(\mathbf{v},z)$ as the contribution from the beam of electrons with velocity $\mathbf{v}$. Using the expressions of $g^{(1/2)}_{B,emg/ref}(\mathbf{v})$ given by Eqs. (11) - (III.3), we can rewrite it

[TABLE]

where, with $s=\pm 1$,

[TABLE]

In the limit $d\rightarrow\infty$, all the exponentials in $L_{emg/ref}$ vanish and we would recover the result for semi-infinite metals; $\mathbf{J}_{B,emg/ref}(z)$ could then be written as a sum of that for two semi-infinite metals. As expected, the surfaces of the film are decoupled in this limit. For thin films, Eq. (23) implies that $\mathbf{J}_{B,emg/ref}(\mathbf{v},z)$ mainly runs along the surface for symmetric modes while normal to it for anti-symmetric modes.

The divergence of $\mathbf{J}_{B}(z)$ can be easily obtained. In the first place we have

[TABLE]

whose Fourier transform is

[TABLE]

with

[TABLE]

Here $\Gamma_{mn}=1+(-1)^{m+n}$, which would vanish identically unless $m$ and $n$ have the same parity. It follows that

[TABLE]

We can write $\mathcal{M}=\mathcal{M}^{+}\bigoplus\mathcal{M}^{-}$, where $\mathcal{M}^{\pm}=\mathcal{M}^{\pm}_{emg}+p\mathcal{M}^{\pm}_{ref}$ operates on the space of $\rho_{\pm}(z)$, with

[TABLE]

In Appendix C, we show that $\mathcal{M}^{\pm}$ is of the order of $kv_{F}/\omega_{p}$.

III.5 Equation of Motion and SPW Solutions

Symmetric and anti-symmetric modes. We proceed to transform Eq. (3) into the equation of motion for $\rho(z)$. In the first place let us show that $\rho_{+}(z)$ and $\rho_{-}(z)$ are strictly decoupled. As is clear from preceding subsections, $\nabla\cdot\mathbf{J}(z)$ and hence the entire left hand side of Eq. (3) are block diagonal with respect to the subspaces respectively spanned by $\rho_{+}(z)$ and $\rho_{-}(z)$. We can prove that $S(z)$ disconnects the subspaces as well. To this end, we Fourier transform $S(z)$ in Eq. (4) to obtain

[TABLE]

Linearly depending on $\rho(z)=\rho_{+}(z)+\rho_{-}(z)$, $J_{z}(z)$ can be split as $J_{z}(z)=J^{+}_{z}(z)+J^{-}_{z}(z)$, where $J^{+/-}_{z}(z)$ denotes the contributions from $\rho_{+/-}(z)$. From their expressions given in preceding sessions, we easily deduce that

[TABLE]

by which we rewrite

[TABLE]

This equation allows us to organize $S_{m}$ in the form of a column vector $\mathcal{S}=\mathcal{S}_{+}\bigoplus\mathcal{S}_{-}$, where $\mathcal{S}_{+,l}=S_{2l}$ contains all the elements $m=2l$ with $l=0,1,...$, while $\mathcal{S}_{-,l}=S_{2l+1}$ contains all the elements $m=2l+1$. As such, the symmetric and anti-symmetric modes belong to different sectors and are strictly decoupled. We can write

[TABLE]

where $\mathbb{E}_{+,l}=1-\delta_{l,0}/2$ and $\mathbb{E}_{-,l}=1$.

Equation of motion. The equation of motion is obtained by Fourier transforming Eq. (3) and using Eqs. (20), (27) and (30) as well as (33). We find

[TABLE]

where the matrix reads

[TABLE]

Here the column vectors are defined by $\rho^{+}_{l}=\rho_{2l}$ and $\rho^{-}_{l}=\rho_{2l+1}$. We can rewrite

[TABLE]

where $\mathcal{G}^{+/-}=\mathcal{G}^{+/-}_{D}+\mathcal{G}^{+/-}_{B}$ is a row vector. We have

[TABLE]

where $G^{+/-}(k,q;\bar{\omega})=G^{+/-}_{D}(k)+G^{+/-}_{B}(k,q;\bar{\omega})$, with

[TABLE]

which is comparable to the counterpart for semi-infinite metals, and

[TABLE]

where $L(v_{x},v_{z},k,q,\bar{\omega},s)=L_{emg}(v_{x},v_{z},k,q,\bar{\omega},s)+p~{}L_{ref}(v_{x},v_{z},k,q,\bar{\omega},s)$, with $L_{emg/ref}$ given by Eqs. (24) and (25).

SPWs as localized solutions. Two types of solutions exist to Eq. (34), depending on whether $J^{+/-}_{z}(0)$ vanishes or not. SPWs are described by solutions with $J^{+/-}_{z}(0)\neq 0$. These solutions represent localized surface waves, for which the equation can be directly solved. We obtain

[TABLE]

which involves no approximations.

Let us write the solution as $\bar{\omega}=\omega_{s}+i\gamma_{0}$ and hence the SPW eigen-frequency is given by $\omega=\omega_{s}+i\gamma$ with $\gamma=\gamma_{0}-1/\tau$. One can show that $\omega_{s}+i\gamma_{0}$ always occurs with $-\omega_{s}+i\gamma_{0}$, in accord with the fact that $\rho(\mathbf{x},t)$ is real-valued. We shall take $\omega_{s}\geq 0$ for definiteness.

Dropping $\mathcal{M}^{+/-}$ as an approximation, the equation becomes

[TABLE]

In addition, we have

[TABLE]

Notably, $\tau$ is not explicitly involved in any of the above equations, implying that the value of $\bar{\omega}$ does not depend on $\tau$.

III.6 Approximate and Numerical Solutions

Hydrodynamic/Drude limits. The hydrodynamic model is attained when the surface-ballistic effects, synthesized in the quantity $G^{+/-}_{B}(k,q;\bar{\omega})$, are ignored in total and the bulk plasma wave dispersion is taken as given by Eq. (21), i.e. $\Omega(k,q;\bar{\omega})\approx\omega^{2}_{p}+(3/5)\left(k^{2}+q^{2}\right)v^{2}_{F}$. In the Drude model, the dispersion is also ignored. In both models, $\bar{\omega}$ is real-valued and Im$(\omega)=-1/\tau$. Solving Eq. (39) without $G^{+/-}_{B}(k,q;\bar{\omega})$, for large $d$ we obtain $\omega^{+/-}_{s}=\omega^{+/-}_{s0}$, with $\omega^{+/-}_{s0}=\left(\omega_{p}/\sqrt{2}\right)\sqrt{1\pm e^{-kd}}$ for the symmetric/antisymmetric modes of SPWs. Note that the bulk wave frequency always lies above the SPW frequency and hence the factor $1/(\Omega^{2}(k,q;\omega_{p})-\bar{\omega}^{2})$ never develops a pole near $\omega_{s}$: SPWs can not decay via bulk waves.

Approximate solutions. We can solve (39) approximately. To the lowest order in $\gamma_{0}/\omega_{s}$, we may determine $\omega_{s}$ by approximating the real part of (39) as follows

[TABLE]

The as-obtained $\omega_{s}$ is then substituted in the imaginary part of Eq. (39) to get $\gamma_{0}$. We find

[TABLE]

which can be brought into a rather simple form if we take $\Omega(k,q;\omega_{s})\approx\omega_{p}$ and $\omega^{2}_{s}/\omega^{2}_{p}\sim 1/2$ for $kd\gg 1$. We get

[TABLE]

with $J^{+/-}_{z}(0)$ evaluated by Eq. (35) with $\omega_{s}$ in place of $\bar{\omega}$. This relation can also be established by an energy analysis, see Sec. IV. By virtue of the relation that $\mbox{Im}\left[G^{+/-}(k,q;\omega_{s})\right]+\mbox{Im}\left[G^{+/-}(k,q;-\omega_{s})\right]=0,$ the same Im$(\bar{\omega})$ exists for $-\omega_{s}$, as anticipated from the fact that charge density waves are real-valued waves.

To make progress, we need to evaluate $G^{+/-}_{B}(k,q;\bar{\omega})$. Writing the integration in Eq. (37) in spherical coordinates and performing it over the magnitude of $\mathbf{v}$, we arrive at

[TABLE]

where we have written $\mathbf{v}=v(\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta)$ and $\tilde{L}(v_{x},v_{z},k,q,\bar{\omega},s)=\left(1\mp e^{i\tilde{\omega}d/v_{z}}\right)~{}L(v_{x},v_{z},k,q,\bar{\omega},s).$ We expand all factors other than $e^{i\tilde{\omega}d/v_{z}}$ in $L(v_{x},v_{z},k,q,\bar{\omega},s)$ into a series of $kv_{F}/\bar{\omega}$ and retain only the leading term. We find

[TABLE]

Upon being formally integrated over $\varphi$, the integral in Eq. (44) ends up in this form, $\int^{1}_{0}dr~{}\left(L_{0}(r)+e^{ir_{0}/r}L_{1}(r)+e^{2ir_{0}/r}L_{2}(r)\right),$ where only the dependence on $r=\cos\theta$ is explicitly noted down in the integrand and $r_{0}=\bar{\omega}d/v_{F}\gg 1$. As $e^{ir_{0}/r}$ and $e^{2ir_{0}/r}$ are rapidly oscillating functions whereas $L_{1,2,3}(r)$ are slowly varying functions, $L_{0}(r)$ is the dominant contribution to the integral. We neglect other contributions and obtain

[TABLE]

This expression explicitly shows that Im$\left[G^{+}_{B}(k,q;\omega_{s})\right]=~{}$Im$\left[G^{-}_{B}(k,q;\omega_{s})\right]<0$, leading to $\gamma_{0}>0$ by virtue of Eq. (42).

It follows that $\mbox{Re}\left[G^{+/-}(k,q;\omega_{s})\right]=(2/d)(\omega^{2}_{p}/4\pi)~{}k~{}\left(1\mp e^{-kd}\right)~{}(1+p)/2.$ Substituting this in Eq. (41) and converting the sum therein into an integral for large $d$, we get $\omega^{+/-}_{s}/\omega_{p}=\sqrt{1-\left(1\pm e^{-kd}\right)(1+p)/4}.$ See that $\omega^{+/-}_{s}$ depends on surface properties via the parameter $p$. Only for $p=1$ would the conventional value, $\omega_{p}/\sqrt{2}$, be recovered. For $p=0$ and at large $kd$, $\omega_{s}=\left(\sqrt{3}/2\right)\omega_{p}$ is slightly larger than the former. It is notable that, $\omega^{+}_{s}$ remains finite even for $k=0$, in distinct contrast with the Drude model. The reason is simple: in Drude model no electric field could exist in the metal for symmetric modes at $k=0$, while in our theory, due to a spatial spread of charge density, the electric field does not vanish. The same conclusion applies to semi-infinite metals.

To estimate $\gamma_{0}$ by Eq. (42), we take in Eq. (46) $\int^{1}_{0}dr~{}\frac{q^{2}v^{2}_{F}r^{3}}{\bar{\omega}^{2}-q^{2}v^{2}_{F}r^{2}}\approx\left(qv_{F}/2\bar{\omega}\right)^{2}$ for simplicity. Thus, $\mbox{Im}\left[G^{+/-}_{B}(k,q;\omega_{s})\right]\approx-(2/d)(\omega^{2}_{p}/4\pi)(3/4)(v_{F}/\omega_{s})~{}(1+p)~{}q^{2},$ which is then plugged in Eq. (42) to produce $\gamma_{0}\sim(3/2\pi)(\omega_{p}/\sqrt{2})/(1\mp e^{-kd})\approx 0.35~{}\omega_{p}.$ In obtaining this expression, we have put $\sum_{l}\frac{\left(q^{+/-}_{l}\right)^{2}}{\left(q^{+/-}_{l}\right)^{2}+k^{2}}\approx q_{c}d/2\pi$ with $q_{c}v_{F}\sim\omega_{p}/\sqrt{2}$. Landau damping has been excluded here, as the approximation only takes the real part of $\Omega(k,q;\omega_{s})$.

Numerical solutions. We can also accurately solve Eq. (39) numerically. The results are displayed in Fig. 1 (a), (b) and (c). A comparison with the approximate solution is not direct, because the approximate solution has excluded while the numerical solution has automatically taken care of Landau damping. It is stressed that, the numerical solutions do not depend on the value of $q_{c}$, provided it is large enough – in excess of $k_{s}$.

IV Energy Conversion with surfaces

In this section, we show that the surface plays a critical role in the energy conversion of bounded systems. While it might be straightforward to handle this issue if the surface potential $\phi_{s}(\mathbf{x},t)$ is exactly known, it is less clear otherwise. Here we derive from Eq. (1) a generic equation that governs the evolution of the electrostatic potential energy, denoted by

[TABLE]

of the system, dispensing with the need to know $\phi_{s}(\mathbf{x},t)$. We then use it to furnish another proof of Eq. (43). For this purpose, we multiply Eq. (1) by $\phi(\mathbf{x},t)$ and integrate it over space to obtain

[TABLE]

where $P^{(1)}(t)=\int d^{3}\mathbf{x}~{}\mathbf{J}(\mathbf{x},t)\cdot\mathbf{E}(\mathbf{x},t)$ is no more than the work done by the electric field on the electrons per unit time and

[TABLE]

It is evident that $P^{(2)}(t)$ signifies the work done by the surface on the electrons per unit time: electrons impinging toward the surface may lose their momentum. As far as we are concerned, this term and its consequences have hitherto not been discussed in existing work. We can translate Eq. (47) into the following, see Appendix A for details,

[TABLE]

where the integral is extended over the metal. If the phase of $\rho(z)$ is global, i.e. independent of $z$, Eq. (49) holds valid even without the abbreviated term.

Now we show how Eq. (43) can also be reached from Eq. (49). Neglecting Landau damping, by Eq. (40) we can show that $\rho(z)$ has a global phase. We can then ignore in this equation the terms abbreviated as Re $\rightarrow$ Im without affecting the results. To the zeroth order in $\gamma_{0}$, it is obvious that Re$\left[\mathbf{J}_{D}(z)\right]\cdot\mbox{Re}\left[\mathbf{E}(z)\right]=0$ and $\mbox{Re}\left[\mathcal{J}_{D}(z)\right]\mbox{Re}\left[E_{z}(z)\right]=0$, i.e. diffusive currents do not bear net work from the electric field. As for the surface-ballistic currents, note that $\mathbf{J}_{B}(\mathbf{v},z)$ contains the rapidly oscillating factor $e^{i\tilde{\omega}z/v_{z}}$, which suppresses the term $\int dz~{}\mbox{Re}\left[\mathbf{J}_{B}(z)\right]\cdot\mbox{Re}\left[\mathbf{E}(z)\right]$ by the factor $kv_{F}/\omega_{p}$, echoing the fact that $\mathcal{M}^{+/-}$ can be neglected in Eq. (34). As such, we have

[TABLE]

where in the last equality, we have used the fact that, for either symmetric or anti-symmetric modes $J_{z}(d)\phi(d)+J_{z}(0)\phi(0)=0$, and that Re$\left[J_{z}(0)\right]=\mbox{Re}\left[J_{B,z}(0)\right]$. To evaluate the denominator, we utilize the equation of motion in real space. It can be easily obtained from Eq. (3) with $\nabla\cdot\mathbf{J}(z)\approx\nabla\cdot\mathbf{J}_{D}(z)\approx(i/\bar{\omega})\omega^{2}_{p}\rho(z)$. We find

[TABLE]

As a result, $\int dz~{}\mbox{Re}\left[\rho(z)\right]\mbox{Re}\left[\phi(z)\right]\approx-\frac{\omega_{s}}{\omega^{2}_{p}-\omega^{2}_{s}}~{}\mbox{Im}\left[J_{z}(0)\right]\mbox{Re}\left[\phi(0)\right].$ By substitution, we immediately recover Eq. (43).

V discussions and conclusions

Thus, on the basis of Boltzmann’s equation, we have established a rigorous theory for SPWs in metal films with arbitrary electronic collision rate $1/\tau$. As a key consequence of the theory, we find that there exists a self-amplification channel for SPWs, which would cause the latter to spontaneously amplify at a rate $\gamma_{0}$ if not for electronic collisions. Surprisingly, the value of $\gamma_{0}$ turns out to be independent of $\tau$. The presence of this channel is guaranteed by the causality principle. Whether the system could actually amplify or not depends on the competition between $\gamma_{0}$ and $1/\tau$. If $\gamma_{0}>1/\tau$, SPWs will amplify and the system will become unstable. In our theory, the non-equilibrium deviation $g(\mathbf{v},z)$ refers to the Fermi-Dirac distribution $f_{0}(\varepsilon)$; as such, the instability is one of the Fermi sea. Needless to say, the instability will be terminated once the system deviates far enough from the Fermi sea and settles in a stable state. Clarifying the nature of the destination state is a subject of crucial importance for future study.

One central feature of our theory is the classification of current densities into a diffusive component $\mathbf{J}_{D}(z)$ and a surface-ballistic component $\mathbf{J}_{B}(z)$. This classification is not based on the value of $\tau$ but according to whether the component obeys the (generalized) Ohm’s law or not. Apart from this, these components are also discriminated in other ways. Firstly, they are controlled by different length scales. As it largely follows the local electric field $\mathbf{E}(z)$, the characteristic length associated with $\mathbf{J}_{D}(z)$ is $k^{-1}$. On the other hand, the length for $\mathbf{J}_{B}(z)$ is $v_{F}/\gamma_{0}$, because of simple $z$-dependence. Secondly, they are oriented disparately. $\mathbf{J}_{D}(z)$ is largely oriented normal to $\mathbf{E}(z)$ locally whereas $\mathbf{J}_{B}(z)$ normal to the surfaces – especially for $p$ close to unity. Considering energy conversion, this explains why $\mathbf{J}_{D}(z)$ does not destabilize the Femi sea but $\mathbf{J}_{B}(z)$ does. Thirdly, $\mathbf{J}_{D}(z)$ is a bulk property and exists regardless of the surface; On the contrary, $\mathbf{J}_{B}(z)$ reflects true surface effects and it would disappear without surfaces.

Although our theory applies at finite temperature, our calculation of $\bar{\omega}$ is done only at zero temperature, i.e. we have taken $f_{0}(\varepsilon)$ to be a step distribution. Clarifying the temperature dependence of Im$(\bar{\omega})$ is important for experimental studies of the present theory, because the net amplification/damping rate $\gamma=\gamma_{0}-1/\tau$ can be directly measured. Arguably, $\gamma_{0}$ could bear a different temperature dependence than $1/\tau$. In sufficiently pure samples, in which the residual resistivity is small enough, there might exist a critical temperature $T^{*}$, above which $\gamma<0$ while below it $\gamma>0$. In other words, $T^{*}$ marks the transition of the system from the Fermi sea to a more stable state.

Another problem that needs to be addressed in the future for experimental studies is concerned with the effects of inter-band transitions. In the most experimented materials, such as silver and gold, these transitions are known to have dramatic effects. They not only open a loss channel due to inter-band absorption, but also significantly shift the SPW frequency. Including them in our formalism consists of a simple generalization: in addition to $\mathbf{J}_{D}(z)$ and $\mathbf{J}_{B}(z)$, the total current density $\mathbf{J}(z)$ must now also have a component $\mathbf{J}_{int}(z)$ accounting for inter-band transitions. The equation of motion is obtained by substituting $\mathbf{J}(z)$ in Eq. (3). One may write $J_{int,\mu}(z)=\sum_{\nu}\int dz^{\prime}~{}\sigma_{\mu\nu}(z,z^{\prime};\omega)E_{\nu}(z^{\prime})$, where $\mu,\nu=x,y,z$ and the inter-band conductivity $\sigma_{\mu\nu}$ can in principle be calculated using Greenwood-Kubo formula. In practice, calculating $\sigma_{\mu\nu}$ could be a formidable task even for the imaginably simplest surfaces. Nevertheless, one may argue that $\mathbf{J}_{int}(z)$ primarily affects the properties of bulk waves, namely, $\Omega(k,q;\bar{\omega})$. The causality principle should still protect the amplification channel, though the value of $\bar{\omega}$ may depend on $\tau$. A systematic analysis will be presented elsewhere.

To conclude, we have presented a theory for SPWs in metal films taking into account the unique interplay between ballistic electronic motions and boundary effects, from which it emerges a universal self-amplification channel for these waves. It is expected that the study will bear far-reaching practical and fundamental consequences, which are to be explored in the future. We hope that the work could stimulate more effort on this subject.

acknowledgement

The author is grateful to K. Wakabayashi for some help with numerical computation. He also thanks E. Mariani for enormous support.

Appendix A More about Eqs. (47) - (49)

The not-so-obvious step in proving Eq. (49) is to show that

[TABLE]

For this purpose, we write $\phi(\mathbf{x},t)=\mbox{Re}\left[e^{-i\omega t}\phi(\mathbf{x})\right]$ and $\rho(\mathbf{x},t)=\mbox{Re}\left[e^{-i\omega t}\rho(\mathbf{x})\right]$, where $\phi(\mathbf{x})=e^{ikx}\phi(z)$ and $\rho(\mathbf{x})=e^{ikx}\rho(z)$. Moreover, we put $\omega=\omega_{s}+i\gamma$, $\phi(\mathbf{x})=\phi^{\prime}(\mathbf{x})+i\phi^{\prime\prime}(\mathbf{x})$ and similarly for other complex quantities. By substitution, we find

[TABLE]

However, $\int d^{3}\mathbf{x}~{}\phi^{\prime}(\mathbf{x})\rho"(\mathbf{x})=\int d^{3}\mathbf{x}~{}\phi^{\prime\prime}(\mathbf{x})\rho^{\prime}(\mathbf{x})=0$. Actually, we have

[TABLE]

thus completing the proof.

Let us suppose $\rho(z)$ has a global phase, i.e. $\rho(z)=c~{}n(z)$, where $c$ is a complex constant and $n(z)$ is real-valued. One can show that Eq. (49) can be turned into an equation that involves only $n(z)$, wherein $c$ plays no role. In other words, Eq. (49) can be evaluated by simply pretending $c$ (and $\rho(z)$) to be real. The proof is evident considering the linear relations between $\rho(z)$ and $\mathbf{J}(z)$ and that between $\rho(z)$ and $\phi(z)$ as well as that between $\rho(z)$ and $\mathbf{E}(z)$.

Appendix B Electronic distribution functions

The general solution to Eq. (7) is given by

[TABLE]

where $C(\mathbf{v})$ is an arbitrary integration constant to be determined by boundary conditions. Let $p_{1}$ and $p_{2}$ be the Fuchs parameters for the (uniform) surfaces at $z=0$ and $z=d$, respectively. The boundary condition at $z=0$ is taken that $g\left((v_{x},v_{y},v_{z}>0),z=0\right)=p_{1}~{}g\left((v_{x},v_{y},-v_{z}),z=0\right)$ while that at $z=d$ assumes $g\left((v_{x},v_{y},v_{z}<0),z=0\right)=p_{2}~{}g\left((v_{x},v_{y},-v_{z}),z=0\right)$, both evaluated at $E_{z}(z)=0$. After some algebra, one finds

[TABLE]

The electronic distribution functions presented in the main text in Sec. III are obtained by approximating $\left(1-p_{1}p_{2}e^{i\frac{2d\tilde{\omega}}{v_{z}}}\right)^{-1}\approx 1$ for $v_{z}\geq 0$ and $\left(1-p_{1}p_{2}e^{-i\frac{2d\tilde{\omega}}{v_{z}}}\right)^{-1}\approx 1$ for $v_{z}<0$ in this equation.

Causality principle. It should be pointed out that, in applying the boundary conditions, we have implicitly assumed Im$(\tilde{\omega})\geq 0$; otherwise, we would find unphysical solutions that violate the principle of causality, which states that the number of out-going electrons is determined by the number of in-coming electrons, not otherwise. It is easy to show that, had we assumed Im$(\tilde{\omega})<0$, we would have found the opposite: the number of reflected electrons would be fixed while the number of incident electrons would go to infinity as $p_{1/2}\rightarrow 0$.

Appendix C The matrix $\mathcal{M}^{+/-}$

In the first place, we show that $\mathcal{M}^{+/-}/\omega^{2}_{p}\propto ikv_{F}/\bar{\omega}+...$, where the ellipsis stands for higher order terms in $kv_{F}/\bar{\omega}$. We take the symmetric modes for illustration, as the reasoning can be replicated for the anti-symmetric modes as well. Writing $\int d^{3}\mathbf{v}~{}\Theta(v_{z})=\int^{2\pi}_{0}d\varphi\int^{\pi/2}_{0}d\theta\sin\theta\int^{\infty}_{0}dv^{2}(v/2)$ and integrating over $v$, we find

[TABLE]

To the lowest order in $kv_{F}/\bar{\omega}$, we only need to retain $L^{(0)}$ in the expansion $L^{sym}=\sum^{\infty}_{m=0}L^{(m)}\left(kv_{F}/\bar{\omega}\right)^{m}$. Thus,

[TABLE]

Substituting this back in (56) and approximating

[TABLE]

we arrive at

[TABLE]

where

[TABLE]

Clearly, we have $\mathcal{M}^{+}/\omega^{2}_{p}\sim kv_{F}/\bar{\omega}$, as stated.

We may proceed further If we take

[TABLE]

from which it follows that $\mathcal{M}^{+}_{l,l^{\prime}}\approx M_{0}=-i\omega^{2}_{p}(1/kd)(kv_{F}/\bar{\omega})$, which is a constant. Therefore, $\mathcal{M}^{+}_{l,l^{\prime}}\approx M_{0}\mathbb{Z}_{l,l^{\prime}}$, where $\mathbb{Z}_{l,l^{\prime}}=1$ constitutes a unity matrix. We write, with $\mathcal{W}^{+}_{l,l^{\prime}}=\delta_{l,l^{\prime}}\Omega(k,q^{+}_{l};\bar{\omega})$,

[TABLE]

where $U$ is a similarity transformation that brings $\mathcal{M}^{+}$ and hence $\mathbb{Z}$ to a diagonal form. We have used a tilde to indicate the transformed matrices, e.g. we write $\tilde{\mathbb{Z}}=U\mathbb{Z}U^{-1}$. See that $\tilde{\mathbb{Z}}$ has only one non-vanishing element, whose value amounts to the dimension $N_{c}$ of the matrix. Let it be the $l_{0}$-th element. Then $\tilde{\mathbb{Z}}_{l,l^{\prime}}=N_{c}\delta_{l,l_{0}}\delta_{l^{\prime},l_{0}}$. Obviously, $N_{c}=q_{c}d/2\pi\sim(\bar{\omega}/kv_{F})(kd/2\pi)$. As such, $M_{0}\sim 1/N_{c}$ and $\tilde{\mathcal{M}}^{+}_{l,l^{\prime}}\sim-i(\omega^{2}_{p}/2\pi)\delta_{l,l_{0}}\delta_{l^{\prime},l_{0}}$. Introducing $\tilde{\mathcal{G}}^{+}=\mathcal{G}^{+}U^{-1}$ and $\tilde{\mathbb{E}}_{+}=U\mathbb{E}_{+}$, we can rewrite the equation of motion for the symmetric modes as

[TABLE]

Taking $\mathcal{W}^{+}\approx\omega_{p}\mathbb{I}$ and hence $\tilde{\mathcal{W}}^{+}\approx\omega_{p}\mathbb{I}$, this equation becomes

[TABLE]

The term in the square bracket makes only a contribution of the order of $\sim 1/N_{c}$ and can be neglected for large $N_{c}$.