Tutorial: The equations of electron emission and their evaluation

Jensen, Kevin L.

doi:10.1063/5.0196481

Electron emission and transport through and over potential barriers is an essential process requiring modeling and simulation to meet the design needs and characterization of an exceedingly broad range of technologically important devices and processes. The simulation and description of thermal, field, and photoemission, and the related concerns of space–charge affected electron flow, often make use of specialized formulations developed in the early days of quantum mechanics. Advancements in the utilization of electron sources and particularly the simulation of devices and applications using advanced particle-in-cell and trajectory methods for beam optics codes create a strong need for a pedagogical account of the emission models to ensure correct numerical evaluation of their equations. This Tutorial starts from simple phenomenological accounts and progressively builds to comprehensive models emphasizing straightforward and often rapid calculation. It recommends formulations to supplant the canonical Richardson–Laue–Dushman (thermal), Fowler–Nordheim (field), Fowler–DuBridge (photo), and Baroody (secondary) equations and provides a useful formulation of space–charge affected flow commonly described by the Child–Langmuir relation that takes into account cathode dependence on surface field.

I. INTRODUCTION

The rapidity with which physicists overcame their skepticism about the existence of an elementary particle of negative charge¹ is a measure of how thoroughly an understanding of its properties became a central concern to modern physics and how that understanding unleashed a wave of revolutionary technological innovations beginning with radar that are still ongoing. Current flow through a vacuum was observed and reported upon in 1884, and its origins were a profound mystery: the current appeared to originate from the hot filament of an early incandescent lightbulb, and because no atoms appeared to be involved, early physicists and electrical engineers struggled for an explanation of the “Edison effect.”² In short order, physicists realized that more than just heat could cause current flow from a cathode to an anode: the application of large potential differences near sharp metal points led to field emission in 1897;³ short wavelength light on metal led to photoemission in 1887;⁴ and the bombardment of materials by high energy primary currents could lead to secondary currents, or secondary emission, in 1902⁵ only a few years after Thomson confirmed in 1897 that cathode rays were due to a new kind of particle dubbed “electrons.”⁶ Before Thomson’s monumental confirmation of the particle basis of cathode rays, the nature of those rays was hotly debated, with a contending view holding that they were an ethereal process in the medium in which light was thought to propagate, but their deflection by magnetic fields was like nothing seen before. Instead, Perrin confirmed the material nature of electrons by deflecting them from entering a hole using a magnet, and Thomson went on to confirm their particle nature, evaluate their charge to mass ratio, and show that their physical size was much smaller than the putatively smallest “atoms.” Langmuir then demonstrated that the charge of the electrons, like the ions observed by Child, suppressed the emission process and so limited the amount of current that could be drawn across a parallel plate anode–cathode (AK) gap.⁷

The canonical equations of electron emission for the five aforementioned processes followed rapidly and are now part of the lexicon of electron emission physics,⁸ with the equations generally known by their developers:

thermal emission (Richardson–Laue–Dushman);
field emission (Fowler–Nordheim);
photoemission (Fowler–Dubridge);
secondary emission (Baroody); and
space charge (Child–Langmuir).

The first four were recast using the methods of the new quantum mechanics in conjunction with Sommerfeld’s model of metals.^9,10 The equations are often referred to by the acronyms of the author names (RLD, FN, FD, CL) except for “secondary emission yield” (SEY), the equation of which is referred to as the “Baroody equation” although it is upgraded by the refinement of others. The unification of a subset of the first three provided an elegant means of demonstrating a common basis for understanding of how electrons surmount or tunnel through the barriers at the surface of a metal that impedes their flow into vacuum. That understanding evolved and led to similar equations governing the transport of electrons between dissimilar materials, especially semiconductor interfaces. Photoemission was amended by Spicer^11,12 through the introduction of a three-step model (TSM)^13–15 that augmented considerations behind the earlier emission equations. Secondary emission likewise made use of a TSM model as well.^16–18

Because the methods of electron emission are so varied, cathodes sometimes are compared on performance metrics that are in fact very different because of how the sources are used. It becomes necessary then to tease apart how the electrons are liberated from the technological methods of how conditions favorable to emission arise for each. It is the purpose of the present treatment to put the most often used equations into a common understanding, extend their range of applicability, and examine how the experimental conditions under which electrons are liberated and then propagate affect the application of the emission equations and their generalization. Methods of calculation are central to that concern: as electron sources and their usages become more ingenious, gaps appear in when the equations are appropriate. Filling those gaps makes the comparative simplicity of the canonical emission equations wistfully remembered. A suite of computational methods that restore some of the lost simplicity is, therefore, a primary goal, as the simulation of modern devices using advanced particle-in-cell (PIC) codes¹⁹ often must account for more than one emission mechanism, but the enormous number of emission sites with local variations in operating conditions plus the requirement to repeatedly evaluate emission over millions of time steps, demands that computational expediency be elevated as a priority. To that end, practical computational and coordinated methods are emphasized.

A. Units

Electron emission occurs over time and length scales billions of times smaller than the macro-scale phenomena that utilize International System of Units “ISU” (aka “MKSA” based on the acronym of meters-kilograms-seconds-amps). It is, therefore, useful to use units that are already atomic scale. The hydrogen atom is the simplest collection of interacting charged particles, and so units based on its physics are preferable. To maintain a straightforward transition back to ISU, prefixed metric units are, therefore, the basis for the preferred dimensional quantities here.

Both energy levels and angular momentum in a quantum mechanical treatment of the hydrogen atom are quantized, and so units natural to it or which have

ℏ

built in are preferred. Conjoining the magnitude of the electron charge, or

q

⁠, with a potential

ϕ

allows

q ϕ \equiv V

to be expressed in (eV), the natural unit to express work function

Φ

in. The ground state kinetic energy of the electron in the Bohr model is the Rydberg energy

R_{y} = (1 / 2) m v^{2} = 13.606

eV, and the angular momentum of the electron is in integer units of Planck’s constant

ℏ

⁠: if the velocity of the electron in the ground state of the hydrogen atom is

v \equiv α c

⁠, where

c

is the velocity of light, then

α = \sqrt{\frac{2 R_{y}}{m c^{2}}} = \sqrt{\frac{27.2113 eV}{510 999 eV}} = \frac{1}{137.036},

(1)

and governs the ratio of the electron kinetic energy to its rest energy. The angular momentum of the ground state electron

m v a_{o} = ℏ

then sets the characteristic unit of length, the Bohr radius

a_{o}

⁠, to be

a_{o} \equiv \frac{ℏ}{m α c} = 5.2918 \times 10^{- 2} nm .

(2)

Finally, the permittivity of free space

ε_{0}

⁠, as used in conventional macroscopic electrodynamics, is more convenient to replace by quantities associated with the potential energy of the ground state electron in the hydrogen atom,

\frac{2 R_{y}}{a_{o}} \equiv 4 Q \Rightarrow Q \equiv \frac{α ℏ c}{2} = 0.359 99 eV nm,

(3)

where the desire to have

Q

defined thusly will become apparent when the image charge potential is encountered in Eq. . Last, a characteristic time is given by how long an electron traveling at its average velocity

α c

takes to circle a hydrogen atom, or

\frac{2 π a_{o}}{α c} = \frac{π ℏ}{R_{y}} = \frac{2.067 83 eV fs}{13.6057 eV} = 0.151 983 fs .

(4)

The natural units for atomic phenomena are then units of (eV) for energy, (nm) for length, (fs) for time, and to take the unit charge to be 1 q: in contrast to MKSA, then, these units are labeled (eV-nm-fs-q)

\to

(enfq). In the present work, all quantities of interest will be in (enfq) units, for which the fundamental constants for electron emission expressed in those units are shown in Table I.

TABLE I.

The fundamental constantsa associated with emission in the units of (eV, fs, nm, q), abbreviated as (efnq). The electron charge is (−q). The unit of temperature remains Kelvin [K].

Symbol	Definition	Value	Unit
−q	Electron charge	−1	q
c	Speed of light	299.793	nm/fs
k_B	Boltzmann’s constant	$1 / 11 604 .5$	eV/K
$ℏ$	Planck’s constant/2π	0.658 212	eV fs
$ℏ$ c	Planck × light speed	197.327	eV nm
α	Fine structure constant	1/137.036	…
m c²	Electron rest energy	510 999	eV
a_o	Bohr Radius $ℏ / m α c$	0.052 917 7	nm
Q	$α ℏ c / 4$	0.359 991	eV nm

Symbol	Definition	Value	Unit
−q	Electron charge	−1	q
c	Speed of light	299.793	nm/fs
k_B	Boltzmann’s constant	$1 / 11 604 .5$	eV/K
$ℏ$	Planck’s constant/2π	0.658 212	eV fs
$ℏ$ c	Planck × light speed	197.327	eV nm
α	Fine structure constant	1/137.036	…
m c²	Electron rest energy	510 999	eV
a_o	Bohr Radius $ℏ / m α c$	0.052 917 7	nm
Q	$α ℏ c / 4$	0.359 991	eV nm

^a

NIST Database https://physics.nist.gov/cuu/Constants/ (accessed 10/07/2023).

Finally, conventions regarding both density and current need to be established. When (enfq) units are being used, it is convenient to render density

ρ (x)

as a number density to facilitate its relation to the same term in Schrödinger’s equation. For ease, consider one dimension, although the arguments convey to three dimensions. Defining the product of the unit charge

q

and the electric potential

φ

to be

V (x)

⁠, then Poisson’s equation transitions to

\frac{d^{2}}{d x^{2}} φ (x) = - \frac{q ρ (x)}{ε_{0}} \Rightarrow \frac{d^{2}}{d x^{2}} V (x) = 16 π Q ρ (x) .

(5)

Hand in hand with the introduction of

V (x)

for potential energy, using a force

F (x) = - q E (x)

rather than an electric field

E (x)

proves convenient for describing phenomena involving the behavior of electrons.

Current is the flow of electrons, for which a convenient conversion is

\frac{1 q}{1 fs} = 6241.5 A .

(6)

The passage of a characteristic charge

q

past a surface in a characteristic time

ℏ / R_{y}

is then

\frac{q c}{a_{o}} = 5665.3 (q / fs) = 0.907 67 A .

(7)

The Alfvén–Lawson current

I_{o}

encountered in the beam envelope (aka rms envelope in the paraxial limit) equation^8,20–22 is related to this characteristic current by

1 / α^{2}

⁠, or

I_{o} \equiv \frac{q c}{α^{2} a_{o}} = 17 045 A .

(8)

and is more commonly encountered as

I_{o} = 4 π ε_{0} m c^{3} / q

in treatments using SI units.

The image charge modified potential barrier that electrons must surmount or overcome to be emitted from a metal naturally sets the form of

Q

⁠. If the metal conduction band minimum is taken as the zero point of energy, then the barrier an electron experiences is given by

V (x) = μ + Φ - F x

in one dimension, where

μ

is the chemical potential (or the Fermi level

E_{F}

at

T = 0

K). When an electron departs a metal surface, the metal remains at zero potential and all field lines from the departing electron are normal to that surface. A mathematically elegant way to meet the boundary conditions is to use the method of images.²³ The image charge exerts a force on the electron as it moves away from the surface, with the distance between the electron and its image being

(2 x)

⁠. As a result, the potential energy

V (x)

is modified to become⁸

V (x) = V_{o} - F x - \frac{Q}{x},

(9)

and is called the Schottky–Nordheim (SN) barrier when

V_{o} = μ + Φ

⁠, alternately, the image charge modified field emission barrier. An example of the SN barrier is shown in Fig. 1. Because of the image charge contribution, the barrier maximum occurs at

x_{o} = \sqrt{Q / F}

and has a height of

V (x_{o}) = μ + Φ - \sqrt{4 Q F}

\equiv μ + Φ - Δ Φ

⁠. Different emission mechanisms have different fields: the Schottky factor for various ones is shown for representative fields in Table II.

FIG. 1.

View large Download slide

Comparison of triangular barrier to image charge barrier for $μ = 7$ eV, $Φ = 4.5$ eV, and $F = 4$ eV/nm for an energy $E = μ - 0.5$ eV, where $μ$ is denoted by the green dashed line. Schottky barrier lowering is governed by $Δ Φ = \sqrt{4 Q F} = 2.4$ eV for which $y (μ) = \sqrt{4 Q F} / Φ = 0.533 327$ ⁠.

TABLE II.

Location of image charge barrier maximum $x_{o} = \sqrt{Q / F}$ and Schottky factor $\sqrt{4 Q F}$ in units of (efnq) for representative $E$ -fields (in common units) associated with various emission mechanisms.

Emission	$E$	F (eV/nm)	ΔΦ (eV)	x_o (nm)
Thermal	20 kV/cm	0.0020	0.0537	13.416
Photo	10 MV/m	0.0100	0.1200	5.9999
Photo	200 MV/m	0.2000	0.5366	1.3416
Field	2 GV/m	2.0000	1.6970	0.4243
Field	8 GV/m	8.0000	3.3941	0.2121

Emission	$E$	F (eV/nm)	ΔΦ (eV)	x_o (nm)
Thermal	20 kV/cm	0.0020	0.0537	13.416
Photo	10 MV/m	0.0100	0.1200	5.9999
Photo	200 MV/m	0.2000	0.5366	1.3416
Field	2 GV/m	2.0000	1.6970	0.4243
Field	8 GV/m	8.0000	3.3941	0.2121

A dimensionless measure of how far the emission barrier has been pulled down by the application of an electric field is given by the ratio of the Schottky factor

Δ Φ

to the height of the barrier above the tunneling electron energy

E

⁠. The convention here, in keeping with the nomenclature of Murphy and Good (MG),²⁴ is to designate that ratio by the symbol

y

⁠, but because this factor shall be generalized below to an electron at an arbitrary energy

E

instead of the energy

μ

corresponding to the Fermi level, then

\begin{matrix} y (E) = \frac{\sqrt{4 Q F}}{μ + Φ - E}, \\ y (μ) = \frac{\sqrt{4 Q F}}{Φ} \end{matrix}

(10)

for the Schottky–Nordheim barrier, for which

y (μ)

corresponds to the

y

of Murphy and Good, and corresponds to the ratio of the Schottky factor to the work function. Observe that if

V_{o}

is defined as the constant part of Eq. (or

μ + Φ

for the SN barrier) and

V (x_{o})

is the maximum of the barrier (or

μ + ϕ

for the SN barrier) and

x = x_{o}

is the location of the maximum, then

y (E) \equiv \frac{V_{o} - V (x_{o})}{V_{o} - E} \equiv \frac{Δ V}{V_{o} - E},

(11)

a form that will be preferable below as it can be generalized to more than just the SN barrier:

y (E)

is the ratio of how much the barrier is lowered through the application of a field compared to how far (in energy) the tunneling electron is below the barrier maximum. In the contrasting but equally applicable formulation of Forbes,²⁵

y {(μ)}^{2} \equiv F / F_{R} \equiv f

(the “scaled field” in Forbes’s nomenclature) holds a special place as the ratio between the applied

F

and a reference

F_{R}

that pulls the SN barrier down to the Fermi level

μ

⁠, and so

y (μ) = 1

represents an upper bound to the applicability of the Fowler–Nordheim equation.²⁶

The evaluation of

μ

can be approximated by the free electron model of Sommerfeld in 1927 that treats the electrons of a metal as a non-interacting gas, the distribution of which follows Fermi–Dirac statistics. A quantum mechanical state can hold either one electron or none (two electrons if spin is included, with spin represented by a factor

g = 2

⁠). The electrons are kept in thermal equilibrium with each other through scattering. At zero temperature, all of the lowest energy states are filled, with the highest filled energy state corresponding to the Fermi energy

E_{F} \equiv ℏ^{2} k_{F}^{2} / 2 m

⁠. At finite temperature,

E_{F}

is replaced by

μ

⁠, and the Fermi–Dirac distribution takes the form

f_{F D} (k) \equiv \frac{g}{1 + \exp [β (E_{k} - μ)]},

(12)

where

β \equiv 1 / k_{B} T

⁠,

T

is temperature measured in Kelvin, and

E_{k} = ℏ^{2} k^{2} / 2 m

⁠.

\vec{k}

is a wave vector such that

ℏ \vec{k}

is the momentum: for an electron in a one-dimensional well of width

L

with infinite walls, then

k = 2 π n / L

⁠, with the integer

n

referring to the energy level.²⁷ Metals typically have Fermi levels of several eV, and so for temperatures encountered for electron sources,

β μ

is large, e.g., for copper with

μ = 7

eV, then

β μ = 63.8112

for

T = 1273

K;

μ (T)

is then to an excellent approximation constant, and integrations that otherwise give rise to factors of

\exp [- β μ]

can replace that term with zero.^8,28,29 The integration of

f_{F D} (k)

over all phase space gives the number

N

of electrons in a box, where

f_{F D}

⁠, being dependent on the energy of the electron, therefore only depends on the magnitude

k = | \vec{k} |

⁠. As a result, although

μ

depends on

T

⁠, the dependence is weak for metals, and so

μ \leftrightarrow ℏ^{2} k_{F}^{2} / 2 m

can be treated as constant given by its zero temperature value to the leading order. If the box is a cube with sides of length

L

⁠, then

μ

can be evaluated in the infinite

β

(zero temperature) limit conveniently, giving

\begin{matrix} N = \frac{1}{{(2 π)}^{3}} \int d^{3} x \int d^{3} k \frac{g}{1 + e^{β (E_{k} - μ)}}, \\ ρ \equiv \frac{N}{L^{3}} \approx \frac{g}{{(2 π)}^{3}} \int_{0}^{k_{F}} 4 π k^{2} d k = \frac{g k_{F}^{3}}{6 π^{2}}, \end{matrix}

(13)

where

E_{F} = μ = ℏ^{2} k_{F}^{2} / 2 m

and

f_{F D} (k)

is

g

for

| \vec{k} | < k_{F}

but vanishes otherwise. The density of the electron gas is then

ρ = N / L^{3}

= k_{F}^{3} / 3 π^{2}

for

g = 2

⁠, as is well known from elementary statistical arguments.³⁰

As an example, consider copper. Naively assuming that each copper atom donates two electrons to the conduction band, then from the atomic weight (⁠

M = 0.063 546

kg/mole) and the mass density (⁠

d_{o} = 8960

kg/m

^{3}

⁠) of copper, it follows

μ \approx \frac{ℏ^{2}}{2 m} {(\frac{3 π^{2} d_{o}}{M})}^{2 / 3} = 7.044 51 eV,

(14)

and is close to the measured value of 7 eV. In general, such agreement among metals is not as fortuitous, but the exercise is sufficient to indicate that Fermi levels of several eV are typical for metals. For that value of

μ

⁠,

k_{F} = 13.5977

(1/nm), and, therefore,

ρ = 8.4912 \times 10^{22}

cm

^{- 3}

⁠, reinforcing that

ρ

is a number density.

B. Basic emission models

If

Δ N

electrons (each of charge

q

⁠) pass a unit area of

Δ A

= Δ x^{2}

in a time

Δ t

⁠, then the current density is the charge flow past that area and is

J = q \frac{Δ N}{Δ x^{2} Δ t} = q (\frac{Δ N}{Δ V}) \frac{Δ x}{Δ t} \to q ρ (x) v_{x},

(15)

where

Δ x^{3} = Δ V

⁠. In the limit of infinitesimal

Δ x \to d x = v_{x} d t

⁠,

ρ (x)

is the number charge density at the position

x

⁠, and the electrons have velocity

v_{x}

in the

\hat{x}

-direction. The sign of

J

indicates the direction of the current flow; for convenience, it is ignored unless needed. Very general arguments can be used to obtain representative order of magnitude estimates of current densities for thermal, field, and photoemission processes. The approximations are necessarily crude so as to introduce the basic notions.

Below, the reader is cautioned that although an effort was made to keep nomenclature/symbols consistent across various sections, the breadth of the subject matter invariably causes symbol choices to clash and their meaning change from one section to the next: to the extent possible, the symbols have been defined in each section in which they appear, but some symbols are reused as needs demand.

1. Simple thermal emission

Thermal emitters using sintered tungsten dispenser cathodes have surface coatings of barium that reduce the work function to

Φ \approx 2

eV. The cathodes are subject to fields comparable to 0.1 kV/cm (⁠

F = 10^{- 5}

eV/nm). The density of electrons is not that of the electron gas in the tungsten metal, but rather a measure of those electrons that are able to surmount the barrier

V_{o} \equiv ℏ^{2} k_{o}^{2} / 2 m

= μ + Φ - \sqrt{4 Q F} \equiv μ + ϕ

where

ϕ = 1.9962

eV. These electrons will be in the thermal tail of the Fermi–Dirac distribution. Crudely, assume that (1/6) of the electrons travel with an average velocity in any given Cartesian direction (a cube has six faces). The proportion of electrons that then surmount the surface barrier

V_{o}

is

(1 / 6)

the total number of electrons with

k > k_{o}

⁠, or³¹

p \approx \frac{1}{6} \frac{\int_{k_{o}}^{\infty} f_{F D} (k) 4 π^{2} k^{2} d k}{\int_{0}^{\infty} f_{F D} (k) 4 π^{2} k^{2} d k} \approx \frac{m k_{o}}{2 β ℏ^{2} k_{F}^{3}} e^{- β ϕ} .

(16)

Approximate the average velocity of an electron in the thermal tail by

(1 / 2) m {⟨ v ⟩}^{2} = k_{B} T / 2

⁠. Assuming

T = 1273

K,

d_{o} = 102 000

kg/m

^{3}

⁠,

M_{W} = 0.183 84

kg/mole, then

k_{o} = 14.2928

(1/nm), and

⟨ v ⟩ = 0.1389

nm/fs. If each tungsten atom donates one electron to the conduction band, then

μ = 9.1862

eV, and

p = 4.1188 \times 10^{- 11}

⁠, giving

J_{T} \approx q p (\frac{d_{o}}{M_{W}}) ⟨ v ⟩ = 5.7951 {A / cm}^{2} .

(17)

A flat thermionic cathode with a radius of (1/4) cm operating at that temperature is then crudely expected to produce a little over 1 A of current.

2. Simple field emission

The largest practical fields generated on the surface of flat cathodes occur in photoinjectors and is limited to at most a few hundred megavolts per meter,³² so that $F_{o} ≲ 0.2$ eV/nm. Because field emission generally requires $F_{o}$ to be in the GV/m range, practical field emitters such as those in Spindt-type field emitter arrays,³³ tungsten wire,^34–36 and carbon fiber^37,38 field emitters, silicon field emitters,³⁹ or carbon nanotubes^40,41 rely on sharpened features that cause fields at the surface to be significantly increased⁴² such that $F = β_{g} F_{o} \sim 4$ –10 eV/nm, where the field enhancement factor $β_{g}$ can be substantial. Note that Forbes et al.⁴² prefer the notation of $γ$ for field enhancement where the surface field $F$ is related to the background field $F_{o}$ ⁠, preferring $β_{g}$ to relate $F$ to a potential difference $V$ ⁠; for continuity here, the present work retains conventions established in Ref. 8, and because $γ$ as a symbol already takes on several meanings later.

A particle viewpoint treats the velocity of an electron as it emerges from the tunneling barrier as zero. It is incompatible with a wave viewpoint that governs attenuation of the wave function within the barrier and with the principle that the velocity and position cannot be exactly known simultaneously. Therefore, the assertion that the velocity of an electron upon emission is zero is unsustainable. Letting current density be represented as

J = q ρ ⟨ v ⟩

and observing that neither the current density

J

nor the number density

ρ

are zero entails that

⟨ v ⟩ > 0

⁠. A wave packet tunneling through a barrier exhibits an exit velocity comparable to its incident velocity in quantum trajectory (Wigner⁴³ and Bohm⁴⁴) models. For purposes of illustration only, assume a representative wave packet has an average velocity of

⟨ v ⟩ \sim ℏ k_{F} / m

(field emitted electrons originate from near the Fermi level). The density is evaluated from the distribution of electrons having tunneled. Let the tunneling probability be based on the Gamow factor

θ (k)

= (4 / 3) κ (k) L

for a triangular barrier of width

L

⁠, where

κ (k) = \sqrt{k_{o}^{2} - k^{2}}

⁠,

μ + ϕ = ℏ^{2} k_{o}^{2} / 2 m

⁠, and

L \sim Φ / F

as will be shown in Sec. II C 2. The tunneling probability is then

D (k) \approx \exp [- θ (k)]

⁠. The fraction

p

that tunnel when

μ = 7

eV,

Φ = 4.5

eV, and

F = 7.5

eV/nm so that

L = 0.6

nm and

ϕ = 1.2137

eV, is then approximately

p \approx \frac{\int_{0}^{k_{F}} e^{- 4 κ (k) L / 3} (k_{F}^{2} - k^{2}) d k}{\int_{0}^{k_{F}} (k_{F}^{2} - k^{2}) d k} \approx 6.997 98 \times 10^{- 9},

(18)

as evaluated numerically, and so with

k_{F} = 13.5546

(1/nm) and

k_{o} = 14.6828

(1/nm), then

J_{F} \approx q p (\frac{k_{F}^{3}}{3 π^{2}}) \frac{ℏ k_{F}}{m} = 2.255 88 \times 10^{8} {A / cm}^{2} .

(19)

A single field emitter with a tip supporting a notional emission area of 10 nm

^{2}

will then produce a current of 22.5588

μ

A; a

100

\times 100

array of identical such emitters would then presumably produce 0.225588 A. That they cannot emit such high currents is a consequence of emission statistics, space charge, emitter heating, and non-uniform conditions over the array.^45–47 Nevertheless, remarkable current densities have been achieved from small arrays.^48–53

3. Simple photoemission

In the simplest model of photoemission from a metal, an electron near the Fermi level absorbs an incident photon of energy

ℏ ω

and subsequently travels in a random direction specified by the polar angle

θ

relative to the surface normal. If the initial energy of the electron is at the Fermi level

μ

⁠, then the total energy after the absorption of a photon becomes

E = μ + ℏ ω

⁠. If the “normal” energy

E_{x} = E \cos^{2} θ = ℏ^{2} k_{x}^{2} / 2 m

exceeds the surface barrier

μ + ϕ

⁠, then the electron is photoemitted. The excited electron need not be exactly at the Fermi level: if

θ = 0

⁠, then the initial energy can be at

ϕ - ℏ ω

below the Fermi level and photoemission can still occur, in keeping with Einstein’s identification of the work function

Φ

being the minimum energy a photon must possess to cause electron ejection⁵⁴ (in the absence of an electric field at the surface). Let

I_{ω}

be the laser intensity (MW/cm

^{2}

⁠) and

[1 - R (ω)] I_{ω}

the amount absorbed where

R (ω)

is the reflectivity. The number of photoexcited electrons

J_{ω} A Δ t / q

⁠, where

Δ t

is the duration of the laser pulse, is then estimated from the number of photons

[1 - R] I_{ω} A Δ t / ℏ ω

absorbed. Quantum efficiency

Q E

is a measure of that ratio, and so, if the duration of the current pulse is equal to that of the laser pulse (which is a good approximation for pulses much longer than the relaxation time for scattering events of excited electrons in metals), then

\begin{aligned} Q E = (\frac{J_{ω} A Δ t}{q}) (\frac{ℏ ω}{I_{ω} A Δ t}) = \frac{ℏ ω}{q} \frac{J_{ω}}{I_{ω}} \equiv [1 - R (ω)] F_{λ} Ω_{p}, \end{aligned}

(20)

where

F_{λ}

measures the loss of electrons due to scattering events as they travel to the surface and

Ω_{p}

is the fraction of a sphere of electrons in energy space that have a normal energy component in excess of the barrier height at the surface. A crude estimate of

Ω_{p}

is to let it be the ratio of the volume of the Fermi–Dirac sphere with

μ ⩾ E ⩾ μ - (ℏ ω - ϕ)

such that the forward energy is in excess of the barrier height compared to the total sphere. To that end, introduce the quantities

Δ = (ℏ ω - ϕ) / μ

and

k_{w} = k_{F} \sqrt{1 - Δ}

⁠, so that

Ω_{p}

is approximated by

\begin{aligned} Ω_{p} & \approx \frac{2 \int_{k_{w}}^{k_{F}} d k_{x} \int_{0}^{\sqrt{k_{F}^{2} - k_{x}^{2}}} 2 π k_{⊥} d k_{⊥}}{\int_{0}^{k_{F}} 4 π k^{2} d k} \\ = 1 - \frac{1}{2} (2 + Δ) \sqrt{1 - Δ} \approx \frac{3}{8} Δ^{2} \end{aligned},

(21)

to leading order in

Δ

⁠. Observe that in the integrals, the differential volume element of the numerator is

d^{3} k = 2 π d k_{x} k_{⊥} d k_{⊥}

(cylindrical coordinates) and the denominator is

d^{3} k = 4 π k^{2} d k

(spherical coordinates). From the definition of

Δ

⁠, then to the leading order

Q E \propto {(ℏ ω - ϕ)}^{2} .

(22)

For metals like copper under laboratory conditions,

μ = 7

⁠,

Φ = 4.5

eV,

F \approx 0.001

eV/nm, and a commonly used laser wavelength of

λ = 266

nm implies

ℏ ω = 4.661 06

eV. Therefore,

Δ = 0.028 429 6

and

Ω_{p} \approx 3.030 90 \times 10^{- 4}

⁠. Measured QEs are lower by about a factor of 20, giving an idea of the magnitude of

(1 - R) F_{λ}

(reflectivity and losses due to scattering) for metals.

4. Simple secondary emission

A beam of high energy electrons (primaries) striking a material will result in the subsequent emission of electrons from the surface, but there are different kinds.⁸ The primaries may collide with electrons in the material and be redirected back to the surface with roughly half their incident energy (back-scattered primaries), or the primaries may undergo repeated scatterings and yet still have enough energy to escape (intermediate-energy primaries), or the electrons with which the primaries scatter have sufficient energy (⁠ $≲ 20$ eV) to be emitted themselves (low-energy secondaries). It is the last group of interest here for $E ≲ 5$ keV.

As the energy $E$ of the primary electron beam increases, the secondary yield $δ$ from solids initially increases, peaks, and then decreases. The maximum yield $δ_{m}$ occurs at $E_{m}$ ⁠. In general, for the metals, the values of the maximum yield typically range between $0.6 < δ_{m} < 1.6$ ⁠, although there is variation depending on whether the surface of the sample is cleaned after insertion into the measurement chamber and even some variation in reported values.⁵⁵ Semiconductors generally have higher SEY’s, with diamond being the highest.^18,56 For purposes of illustration, consider the data synopsized in Figs. 1 and 2 of Baroody¹⁶ and shown in Figs. 2 and 3. Baroody had shown $δ_{m}$ as a function of $Φ$ to show that they scaled together: the present representation conveys the same message but shows the influence of where the element sits on the Periodic Table, with the green lines intersecting elements in the alkali column.

FIG. 2.

Ratio of SEY δ ( E ) to maximum value δ m for six metals identified in the legend. Metal data were digitally extracted from Fig. 1 of Baroody.16 Curves are based on Eqs. (25) and (26), where the choice of λ is comparatively flexible for the ratios.

View large Download slide

Ratio of SEY $δ (E)$ to maximum value $δ_{m}$ for six metals identified in the legend. Metal data were digitally extracted from Fig. 1 of Baroody.¹⁶ Curves are based on Eqs. (25) and (26), where the choice of $λ$ is comparatively flexible for the ratios.

FIG. 3.

Secondary yield δ m and work function Φ for various metals as a function of atomic number, based on data digitally extracted from Fig. 2 of Baroody16 and rounded to the hundredth decimal place. Φ has been scaled by f ( Φ ) = ⌊ δ m ⌋ + ( Φ − ⌊ Φ ⌋ ) ( ⌈ δ m ⌉ − ⌊ δ m ⌋ ) ] / ( ⌈ Φ ⌉ − ⌊ Φ ⌋ ), where ⌊ x ⌋ = min ( x ) and ⌈ x ⌉ = max ( x ). The elements shown, in order of atomic number, are (Mg, Al, K, Fe, Co, Ni, Cu, Rb, Zr, Mo, Pd, Ag, Cd, Sb, Cs, Ba, Ta, W, Au, Th). Green lines show AN of alkali metals.

View large Download slide

Secondary yield $δ_{m}$ and work function $Φ$ for various metals as a function of atomic number, based on data digitally extracted from Fig. 2 of Baroody¹⁶ and rounded to the hundredth decimal place. $Φ$ has been scaled by $f (Φ) = ⌊ δ_{m} ⌋ + (Φ - ⌊ Φ ⌋) (⌈ δ_{m} ⌉ - ⌊ δ_{m} ⌋)] / (⌈ Φ ⌉ - ⌊ Φ ⌋)$ ⁠, where $⌊ x ⌋$ $= min (x)$ and $⌈ x ⌉ = max (x)$ ⁠. The elements shown, in order of atomic number, are (Mg, Al, K, Fe, Co, Ni, Cu, Rb, Zr, Mo, Pd, Ag, Cd, Sb, Cs, Ba, Ta, W, Au, Th). Green lines show AN of alkali metals.

A simple model of Baroody^8,16 predicts that

δ (E) / δ_{m}

as a function of

E / E_{m}

falls on a universal curve as suggested by Fig. 2. As primaries penetrate a material, they lose energy after collisions with other electrons, suggested to be

\frac{d E}{d x} \approx - \frac{A}{E^{n - 1}},

(23)

where

A

is a proportionality constant and Baroody lets

n = 2

(Whiddington’s law). The decline in energy

E (x)

is then

\int_{0}^{x} d x^{'} = \int_{E_{o}}^{E (x)} \frac{d x}{d E} d E \to E (x) = {[E_{o}^{n} - n A x]}^{1 / n},

(24)

where

E_{o}

is the initial energy of the primary electron. The range

R (E_{o})

is when the primary has lost all of its energy (⁠

E (R) = 0

⁠), and so is found to be

R (E_{o}) = E_{o}^{n} / n A

⁠. If the loss of primaries as a function of

x

due to various mechanisms is parametrically represented as

N (x) = N_{o} \exp (- γ x)

⁠, where

N (x)

is the number of primaries at

x

⁠,

N_{o}

is the initial number of primaries, and

γ

is a decay depth factor, then

\begin{aligned} δ (E_{o}) = - B \int_{0}^{R (E_{o})} \frac{d E}{d x} e^{- γ x} d x = B E_{o} G_{n} [γ R (E_{o})], \\ G_{n} (λ) \equiv e^{- λ} \int_{0}^{1} \exp (λ s^{n}) d s, \end{aligned}

(25)

where

λ \equiv γ R (E_{o}) = γ E_{o}^{n} / n A

(here, it is not a wavelength), and

B = D / Δ E

is taken as the ratio of a surface emission probability

D

with the average energy

Δ E

required to free a core electron from one of the atoms making up the secondary emission material.⁵⁷ The factor

e^{- γ x}

is the probability that a secondary electron generated at a depth

x

into the material can reach the surface without losing energy due to scattering which renders it emission-ineligible.⁵⁸ The function

G (λ)

results from using Eq. and defining

A n x = s^{n} - 1

⁠. For

1 ⩽ n ≲ 2

⁠, a reasonably good approximation is

G_{n} (λ) \approx \frac{1 - e^{- λ n}}{λ n} .

(26)

Next, after evaluating

δ (E_{o})

for a range of

E_{o}

⁠, find the

E_{m}

such that

δ (E_{m})

is the largest value, then form the ratios

E_{o} / E_{m}

and

δ (E_{o}) / δ (E_{m})

⁠. It is found that to a very good approximation, curves defined by these ratios are visually independent of the values of

λ

for a given

n

⁠, and that larger

n

cause

δ (E_{o}) / δ (E_{m})

to decline more rapidly with

E_{o}

⁠. The particular choices

(n, λ) = (1.4, 5)

and

(n, λ) = (2, 2)

are shown in Fig. 2, showing that other values of

n

other than 2 appear to match the data better, even though the particular value of

λ

is flexible.

Later studies by Lye and Dekker⁵⁸ taking into account findings by Young⁵⁹ for secondary emission from aluminum led them to modify the approximation to

d E / d x

and subsequently conclude that the universal curve takes the form of⁸

\frac{δ (E_{o})}{δ (E_{m})} = \frac{1}{ξ^{n - 1}} [\frac{1 - e^{- λ ξ^{n}}}{1 - e^{- λ}}],

(27)

for

ξ = E_{o} / E_{m}

⁠, but a different approach shall be taken in Sec. III E when a more complex secondary emission relation is undertaken.

5. Simple space charge

An excess of charge between a cathode and anode (taken to be metals) produces fields at the cathode that can reduce subsequent emission given that all of the simple emission models depend in one way or another on

F

brought on by the electric field at the surface of the cathode. In a one-dimensional problem, differential elements of charge correspond to sheets of charge, where the charge per unit area of the sheet is the product of the number charge density

ρ (x)

and the differential element

d x

⁠, or

ρ (x) d x

⁠. It is known from standard electrostatics using a Gaussian pillbox argument^23,60 that the field to each side of the charge sheet is constant, independent of distance to the sheet, and proportional to the charge per unit area. By superposition, the field downstream from the cathode surface changes incrementally with each additional sheet, or [assuming the charge in the gap is due to electrons with a charge of

(- q)

]

d F = \frac{q^{2}}{K_{s} ε_{0}} ρ (x) d x,

where

K_{s}

is the dielectric constant of the gap region (⁠

K_{s} = 1

for vacuum) and

ε_{0}

is the permittivity of free space. To make the units match those of the present Tutorial, this equation is rewritten as

\frac{d F}{d x} = \frac{d^{2}}{d x^{2}} V (x) = \frac{16 π Q}{K_{s}} ρ (x),

(28)

where

V (x) = q φ (x)

is a potential energy if

φ (x)

is a potential, and the sign of the equation is a consequence of the negative charge of the electron. From Table I,

16 π Q = 2.899 16 \times 10^{- 27} J m

= 18.0951

eV nm. The solution to Eq. depends on whether number density

ρ (x)

is constant (as occurs for uniform doping in a crystal lattice) or whether current density

J (x)

is, as is relevant to vacuum electronic devices and their electron sources. Consider each in turn.

For constant number density

ρ (x) = ρ_{o}

⁠, then for a gap of width

D

⁠, a double integration yields⁸

\frac{V_{a}}{D} = F_{c} + \frac{8 π Q D}{K_{s}} ρ_{o},

(29)

where

F_{c}

is the field at the boundary between metal and the gap (taken to be filled with a doped semiconductor). The “depletion width”

w

is defined by setting

F_{c} = 0

to zero,

ρ_{o} \to N_{d}

(the doping level), and

D \to w

to obtain

V_{a} = \frac{8 π Q w^{2}}{K_{s}} N_{d},

(30)

where

V_{a}

is the barrier preventing electron injection from a metal into a semiconductor thin film. The shape of

V (x)

is parabolic with the origin at the metal–semiconductor interface and the minimum occurs where

x = w

(the depletion width), beyond which

V (x)

is flat. As an example for parameters chosen to mimic diamond single crystal thin film studies,¹⁷ for an interface barrier of

V_{a} = 4.64

eV, a dielectric constant of

K_{s} = 5.7

⁠, and a boron doping of

N_{d} = 4 \times 10^{13}

atoms/cm

^{3}

⁠, then

w = 8.5487 μ

m, a value which drops to

w = 540.67

nm if

N_{d} \to 10^{16}

atoms/cm

^{3}

⁠.

Current density, understood to be the passage of charge past an area in a given time

Δ t = Δ x / v (x)

⁠, is approximated by

J = q v (x) ρ (x)

⁠, where

v (x)

is the electron velocity in vacuum where

K_{s} \to 1

⁠. Solving Poisson’s equation for the space charge problem^8,61–64 is accomplished by writing

d^{2} V / d x^{2} = d F / d x

= (1 / 2) d (F^{2}) / d V

and, from conservation of energy,

(1 / 2) m v {(x)}^{2}

= V_{m} + V (x) \equiv U (x)

⁠, where

V_{m}

is the cathode potential and for convenience set to 0.

V (x)

is the potential energy of the electron in the anode–cathode (AK) gap: an anode at a higher potential

φ (D)

corresponds to a lower potential energy

V (D) = - q φ (D) = - V_{a}

(where

V_{a}

is taken as a positive quantity to enable convenience later). Poisson’s equation is, therefore, rewritten as

d F^{2} = \frac{16 π Q}{q} J {[\frac{2 m}{U}]}^{1 / 2} d U,

(31)

for which the solution is provided by elementary integration to be

F^{2} - F_{m}^{2} = (16 π Q J \sqrt{2 m} / q) \sqrt{U}

⁠, where

F_{m}

is the

F

value at the surface of the cathode. Using

F \to d U / d x

allows for a second integration to be performed, for

x

ranging from 0 to

D

and

U

from 0 to

V_{a}

⁠. For the special case where the field at the surface of the cathode is

F_{m} = 0

⁠, then the current density

J = J_{C L}

is the Child–Langmuir result and is

J_{C L} (V_{a}, D) = \frac{q}{18 π Q \sqrt{2 m}} \frac{V_{a}^{3 / 2}}{D^{2}} .

(32)

Returning Child’s law to its more familiar SI units, Eq. for

V_{a} = q φ_{a}

becomes

J_{C L} (φ_{a}, D) = \frac{4 ε_{0}}{9 D^{2}} \sqrt{\frac{2 q}{m}} φ_{a}^{3 / 2} .

(33)

As an example, for

V_{a} = 1

eV and

D = 1

μ

m,

J_{C L} = 1.456 74

\times 10^{- 8}

q/fs nm

^{2}

= 233.395 A/cm

^{2}

⁠, corresponding to about 15 electrons every picosecond from a square micrometer area.

II. QUANTUM MECHANICAL CURRENT DENSITY

A. Distributions

A collection of particles distributed with respect to both position and velocity introduces a distribution function

f (\vec{r}, \vec{k}; t)

in statistical mechanics,^30,65–67 a well-known example being the time-independent Maxwell–Boltzmann distribution when classical particles are uniformly distributed spatially. The distribution of quantum particles (fermions and bosons) have as their asymptotic limit the MB distribution, and so the MB distribution often functions as an easier means to confirm statistical arguments, and sometimes (as in the case of the RLD equation) used directly. The canonical electron emission equations assume a time-independent one-dimensional potential barriers

V (x)

⁠. Consequently, the conjugate momentum

k_{x}

can be abbreviated to

k

⁠, and rather than a six-dimensional distribution

f (\vec{r}, \vec{k}, t)

⁠, it is sufficient instead to treat a two-dimensional distribution

f (x, k, t)

⁠, delaying the generalization to six dimensions to when it is needed. Moments of the distribution give number density and current density via

\begin{matrix} ρ (x, t) = \frac{1}{2 π} \int_{- \infty}^{\infty} f (x, k, t) d k, \\ J (x, t) = \frac{q}{2 π} \int_{- \infty}^{\infty} \frac{ℏ k}{m} f (x, k, t) d k . \end{matrix}

(34)

Although

ρ (x, t)

is a number density, the inclusion of

q

in the definition of

J (x, t)

renders it a charge current density as is more convenient. In a phase space representation of quantum mechanics, these equations continue to function, but they are come to by means of constructing a non-local phase space description obtained from the density matrix (the Wigner function^8,68–70). If particle number is conserved (⁠

f

is treated like an incompressible fluid), then

0 = \frac{d}{d t} f (x, k, t) = \partial_{t} f + \dot{x} \partial_{x} f + \dot{k} \partial_{k} f,

(35)

where the Newtonian over-dot notation denotes time derivative of

x

and

k

on the right hand side, and the resulting equation is recognized as Boltzmann’s transport equation.^71,72 Its extension to treat quantum particles causes the last term

\dot{k} \partial_{k} f

to be replaced by a non-local function dependent on

V (x \pm y)

away from the phase space location

x

⁠. The local force is

\partial_{x} V (x) = - ℏ \dot{k}

⁠; as a result, if

V (x)

is independent of time, then integrating both sides over all

k

⁠, assuming

f (x, \pm \infty) \to 0

and using Eq. gives the continuity equation

q \partial_{t} ρ (x, t) + \partial_{x} J (x, t) = 0,

(36)

where, by virtue of the way

J (x)

is defined, a

q

is appended to

ρ (x)

⁠.

For classical point particles where

x

and

k

are independently treated, Eq. provides the foundation for treating

J

\sim q ρ v

as in the “simple” equations of Sec. I B. At the atomic scale, however, particles such as electrons have a wavelength

λ

governed by the ratio of Planck’s constant with the particle momentum (the de Broglie relation). Electrons at the Fermi level, therefore, satisfy

λ = 2 π / k_{F}

⁠: for metals, this length is typically sub-nanometer and, therefore, comparable to the width of the barrier for field emission conditions. The nature of transport past (or through) the barrier is, therefore, inherently quantum mechanical and requires solutions to Schrödinger’s equation which, in one-dimension, is

(- \frac{ℏ^{2}}{2 m} \frac{d^{2}}{d x^{2}} + V (x)) ψ_{k} (x) = E_{k} ψ_{k} (x),

(37)

where

E_{k} \equiv ℏ^{2} k^{2} / 2 m

⁠, and

k

is the conjugate variable to position

x

⁠, where

ℏ k

is momentum, and

ψ_{k} (x)

is an energy eigenstate: Transfer Matrix Approach (TMA) solutions to Schrödinger’s equation⁸ have analogs to perhaps more familiar electromagnetic methods.⁷³ A weighted sum over eigenstates gives the number density in the quantum formulation,

ρ (x) = \sum_{n = 1}^{\infty} C_{n} | ψ_{n} (x) |^{2},

(38)

where

C_{n}

is the weighting factor for the

n

th eigenstate characterized by

E_{n} = ℏ^{2} k_{n}^{2} / 2 m

if the eigenstates are discrete. When extended to the continuum,

C_{n}

is replaced by a weighting factor

f (k)

(the notation of which is intentional). In the simple model of a plane wave incident on a barrier, e.g., the conventional one being a rectangular barrier separating two flat (in energy) regions, then the outgoing plane wave has the form

ψ_{k} (x) = t (k) e^{i k x}

⁠, where

t (k)

is the transmission coefficient from which the transmission probability

D (k) = | t (k) |^{2}

is constructed. It, therefore, follows

ρ (x) \equiv \frac{1}{2 π} \int_{- \infty}^{\infty} D (k) f (k) d k,

(39)

a form appealing for its similarity to Eq. . From quantum mechanics,^8,27,74 the probability current density associated with an eigenstate is^8,75

j_{k} = \frac{ℏ}{2 m i} [ψ^{†} \partial_{x} ψ - ψ \partial_{x} ψ^{†}] = \frac{ℏ k}{m},

(40)

and so, weighting the eigenstates and summing over them gives the equation for current density as

J (F, T) = \frac{1}{2 π} \int_{- \infty}^{\infty} \frac{ℏ k}{m} D (k) f (k) d k,

(41)

which is related to the Tsu–Esaki formula⁷⁶ for tunneling when energy

E = ℏ^{2} k^{2} / 2 m

is parabolic in momentum, even though the “derivation” here is not rigorous and merits examination.⁷⁷ Equation is the defining equation, though, for current density that allows the convenient derivation of the canonical emission equations and is used exclusively herein, although some authors^24,78–81 work in energy

E

rather than

k

as here. The velocity term

ℏ k / m

is readily understood, but the weighting factor (or “supply function”)

f (k)

and the transmission probability

D (k)

require specification.

B. Supply function

For energy parabolic in momentum, the factor

(ℏ k / m) f (k) d k

could be re-expressed as

N (E) d E

⁠, allowing for a useful means of evaluating current density in terms of energy integrals.^1,80 Then,

N (E) d E

is referred to as the normal energy distribution giving number of electrons having a normal energy between

E

and

E + d E

⁠, alternately the incident normal energy distribution,⁸¹ that strike a unit area of the surface in unit time. However, the reliance of the present treatment on wave number

k

rather than energy

E

reflects an emphasis on providing connections to the calculation of the quantum mechanical tunneling probability using the conjugate variable to

x

as required (switching to an energy integration is, by comparison, a convenient mathematical trick motivated by the Sommerfeld model of an electron gas) using Airy-Transfer Matrix Approaches (TMAs). Consequently, a naming convention introduced in the early literature^24,79 is kept so that the weighting factor

f (k)

is termed a “supply function.” It is seen that the supply function, as defined herein, is, in fact, the Fermi–Dirac distribution integrated over the transverse momentum

k_{⊥}

components, and so

f (k_{x}) \equiv \frac{g}{{(2 π)}^{2}} \int_{0}^{\infty} \frac{2 π k_{⊥} d k_{⊥}}{1 + \exp [β_{T} (E_{x} + E_{⊥} - μ)]},

(42)

where

E_{x} = ℏ^{2} k_{x}^{2} / 2 m

⁠,

E_{⊥} = ℏ^{2} k_{⊥}^{2} / 2 m

⁠, and the total energy is their sum. The factor

β_{T}

will be referred to as the energy slope factor for temperature. The integral is analytic of the form

\int d u / (1 + e^{u}) = \ln (1 + e^{- u})

⁠, so that with

g \to 2

⁠,

f (k_{x}) = \frac{m}{π β_{T} ℏ^{2}} \ln {1 + \exp [β_{T} (μ - E_{x})]},

(43)

after which, in one dimension, the

x

-subscript is again ignored. The canonical emission equations generally rely on two limits:

thermal limit: when $β_{T} Φ ≫ 1$ and $k > k_{F}$ ⁠, then

$f (k) \to \frac{m}{π β_{T} ℏ^{2}} e^{β_{T} (μ - E)},$
(44)
field limit: when $T \to 0$ (i.e., $β_{T} \to \infty$ ⁠) then $f (k)$ vanishes for $k > k_{F}$ ⁠, but when $k < k_{F}$ ⁠, then

$f (k) \to \frac{1}{2 π} (k_{F}^{2} - k^{2}),$
(45)
which recovers Eq. (13) by

$\frac{1}{2 π} \int_{- k_{F}}^{k_{F}} f (k) d k = \frac{k_{F}^{3}}{3 π^{2}} .$
(46)

C. Transmission probability

Exact methods to calculate the wave function for an electron encountering a barrier exist based on Transfer Matrix Approach (TMA) methods,⁸ but for the purposes of developing the canonical emission equations, a simpler approach suffices. The wave function

ψ (x, t)

is a solution to the time-dependent Schrödinger equation

[- \frac{ℏ^{2}}{2 m} \frac{\partial^{2}}{{\partial x}^{2}} + V (x)] ψ (x, t) = i ℏ \frac{\partial}{\partial t} ψ (x, t) .

(47)

Let the wave function be the product of a real coefficient

R (x, t)

and the exponential of a phase

i S (x, t)

⁠, or

ψ = R e^{i S}

⁠. When inserted into Eq. , then a variety of terms appear, some real and others imaginary. The sum of all the real terms must vanish separately from the sum of all the imaginary terms, and so

\begin{matrix} 0 = \frac{ℏ^{2}}{2 m} [{(\partial_{x} S)}^{2} - R^{- 1} \partial_{x}^{2} R] + V (x) + ℏ \partial_{t} S, \\ 0 = i ℏ R^{- 1} \partial_{t} R + i \frac{ℏ^{2}}{2 m R} [2 \partial_{x} R \partial_{x} S + R \partial_{x}^{2} S] . \end{matrix}

(48)

Make use of a few replacements. The density is

ρ = | ψ (x, t) |^{2} = R^{2}

⁠. From Eq. , current density is

j = ρ (ℏ / m) \partial_{x} S

⁠, a form which identifies

(ℏ / m) \partial_{x} S

as a velocity term. Define the quantum potential⁸² by

Ω = - (ℏ^{2} / 2 m) R^{- 1} \partial_{x}^{2} R

and replace

\partial_{x} S

in favor of

j

⁠,

\begin{matrix} 0 = \frac{i ℏ}{2 ρ} (\partial_{t} ρ + \partial_{x} j), \\ 0 = ℏ \partial_{t} S + \frac{ℏ^{2}}{2 m} {(\partial_{x} S)}^{2} + Ω (x) + V (x), \end{matrix}

(49)

where the first equation is seen to be a quantum mechanical continuity equation analogous to Eq. . Using

i ℏ \partial_{t} ψ = E ψ

to obtain

\partial_{t} S = (i / R) \partial_{t} R

⁠, and restricting attention to steady state conditions for which

\partial_{t} ρ = \partial_{x} j = 0

finally gives the single equation

{(\partial_{x} S)}^{2} = \frac{2 m}{ℏ^{2}} [E - V (x) - Ω (x)] .

(50)

If the density

ρ

is slowly varying, then the quantum potential is negligible. Finally, when the energy

E

is below the barrier

V (x)

(tunneling), then

\partial_{x} S

is imaginary, and so the wave function exponentially decays in passing through the barrier. The measure of decay is given by integrating

\partial_{x} S

between the zeros of

V (x_{\pm}) - E = 0

⁠, where

x_{\pm}

is defined by this equation. If the ratio of the density on the right hand side of the barrier with the density on the left side defines the transmission probability, then that identifies

D (E) = e^{- θ (E)},

where

θ (E) = 2 \frac{\sqrt{2 m}}{ℏ} \int_{x_{1}}^{x_{2}} \sqrt{V (x) - E} d x,

(51)

which is commonly called the Gamow factor of tunneling.⁸³ The limits of integration are defined by the smallest two roots of

V (x_{j}) - E = 0

if more than two exist, or

j \in (1, 2)

if only two roots exist (as for the rectangular, triangular, quadratic, and SN barriers). Compared to conventional methods in electron emission theory, the following modifications are utilized here. First, the Kemble form of the transmission probability^84,85 is relied upon in the master equation, Eq. , for current density

D (k) \to \frac{1}{1 + e^{θ (k)}} .

(52)

For parabolic barriers, the Kemble form is exact, and for general barriers, it is superior to the conventional form of

D (k) \sim e^{- θ}

⁠. Second, a shape factor methodology^8,86 will be exclusively used in evaluating the Gamow factor, now defined as

\begin{matrix} θ (E) \equiv 2 σ [y (E)] κ (E) L (E), \\ κ (E) \equiv \frac{1}{ℏ} \sqrt{2 m [V_{o} - E]} = \sqrt{k_{o}^{2} - k^{2}}, \\ L (E) \equiv x_{2} (E) - x_{1} (E), \end{matrix}

(53)

for which the shape factor is defined by⁸⁶

σ [y (E)] = \int_{x_{1}}^{x_{2}} {[\frac{V (x) - E}{V (x_{o}) - E}]}^{1 / 2} \frac{d x}{L (E)},

(54)

where

y (E)

relies on Eq. for general barriers and

V (x_{o})

is the maximum height of the barrier which occurs at

x = x_{o}

⁠. The shape factors of several common barriers are analytic.

1. Rectangular barrier

The rectangular barrier for which

V (x) = V_{o} Θ [{(L_{o} / 2)}^{2} - x^{2}]

is the simplest, for which

x_{\pm} = \pm L_{o} / 2

⁠, where

L_{o}

is understood to be where

V (x)

intersects the

x

axis and

Θ (x)

is the Heaviside step function. It follows

σ_{□} [y (E)] = \int_{- L_{o} / 2}^{L_{o} / 2} \frac{d x}{L_{o}} = 1.

(55)

As a consequence, the Gamow factor for the rectangular barrier is

θ_{o} \equiv 2 κ (E) L_{o}

⁠, where

κ^{2} = k_{o}^{2} - k^{2}

or, in terms of energy

E = ℏ^{2} k^{2} / 2 m

⁠,

θ_{o} (E) = \frac{2 L_{o}}{ℏ} \sqrt{2 m (V_{o} - E)} .

(56)

2. Triangular barrier

The triangular barrier for which

V (x) = (V_{o} - F x) Θ (x)

was considered by Fowler and Nordheim⁸⁷ as it allows for an exact solution of

D (E)

in terms of Airy functions without relying on the Gamow factor.⁸ Nevertheless,

σ_{▹} [y (E)] = \int_{0}^{L (E)} {[\frac{V_{o} - E - F x}{V_{o} - E}]}^{1 / 2} \frac{d x}{L (E)} = \frac{2}{3},

(57)

where

L (E) = (V_{o} - E) / F

⁠. Consequently, the triangular barrier Gamow factor is

θ_{▹} = (4 / 3) κ (E) L (E)

as found by Fowler and Nordheim using asymptotic forms of Airy functions [compare the argument of the exponential of

D_{F N} (k)

in Eq. from Appendix A].

3. General symmetric barrier

The class of symmetric barriers for

V_{n} (x) = V_{o} [1 - {(\frac{| x |}{L_{o}})}^{2 / n}],

(58)

have constant shape factors.⁸⁶ Introducing the factor

λ_{o}

= V_{o} / (V_{o} - E) = k_{o}^{2} / κ {(E)}^{2}

such that

x_{2} = L_{o} λ^{- n / 2}

⁠, then

\begin{aligned} σ_{n} & = \frac{1}{L (E)} \int_{0}^{L (E)} {[1 - λ_{o} {(\frac{x}{L (E)})}^{2 / n}]}^{1 / 2} d x \\ = \int_{0}^{1} {(1 - s^{2 / n})}^{1 / 2} d s = \frac{Γ (3 / 2) Γ [(n / 2) + 1]}{Γ [(n + 3) / 2]}, \end{aligned}

(59)

for which

σ_{0} = 1 = σ_{□}

is the rectangular barrier,

σ_{2} = 2 / 3

for the isosceles triangle is equivalent to

σ_{▹}

(a right triangle necessarily has the same shape factor as an isosceles triangle formed by reflecting the right triangle about its vertical side), the parabolic barrier behind the Kemble approximation for which

σ_{1} \equiv σ_{\cap} = π / 4

⁠, and

σ_{3} = 3 π / 16

⁠. As

n

increases,

σ_{n}

monotonically decreases⁸⁸ crudely as

1 / \sqrt{n}

⁠.

4. Trapezoidal barrier

Metal–insulator–metal (MIM) tunnel junctions,⁸⁹ and the related metal–oxide–semiconductor (MOS) structures are a trapezoidal barrier that is an extension of the triangular barrier of Sec. II C 2 but which allows for the first of the

y

-dependent shape factors to be considered. As dimensions of nano-electronic devices continue to shrink, MIM barriers are increasingly important in nano- and subnano-scale junctions and contacts.^90,91 Image charge effects complicate the MIM model beyond the Schottky–Nordheim (SN) complications^86,89 considered in Sec. II C 5. The barrier is represented as

V_{⋄} (x) \equiv (V_{o} - F x) Θ (x) Θ (L_{o} - x),

(60)

where

L_{o}

is the fixed width of the barrier: consequently, the energy dependence of

σ_{⋄}

is determined by the behavior of

κ (E)

⁠. Because

y (E)

defined by Eq. is always zero, introduce the dimensionless factor

λ (E)

defined by

λ (E) = \frac{L_{o} F}{V_{o} - E},

(61)

in terms of which

σ_{⋄}

for

V_{o} > E

is

σ_{⋄} (λ) = \frac{2}{3 λ} [1 - {(1 - λ)}^{3 / 2}],

(62)

and shown in Fig. 4. When

E > V_{o}

or

λ > 1

⁠, then

σ_{⋄} \to σ_{▹}

because

L_{o} \to L (E)

from Sec. II C 2.

5. Schottky–Nordheim barrier

The Schottky–Nordheim (SN) barrier is the basis for the Richardson–Laue–Dushman and Fowler–Nordheim equations as most commonly used: it, like the metal–insulator–metal (MIM), trapezoidal, or depletion barriers,⁸⁶ does not have a constant shape factor, and so requires greater attention. For the SN barrier [Eq. ],

V (x) - E = \frac{- Q + (V_{o} - E) x - F x^{2}}{x} .

(63)

FIG. 4.

View large Download slide

Comparison of trapezoidal barrier $σ_{⋄} (λ)$ [Eq. (62)] to triangular $σ_{▹}$ $= 2 / 3$ and rectangular barrier $σ_{□} = 1$ ⁠.

The numerator is a quadratic polynomial for which two real roots exist corresponding to

V (x_{j}) - E = 0

between which exists a maximum at

x_{o} = \sqrt{Q / F}

⁠: identifying the first root as

x_{1}

and the second as

x_{2}

⁠, letting

Δ (E) \equiv V_{o} - E

⁠, and using

y (E) = \sqrt{4 Q F} / Δ

for the SN barrier results in (where now

j

is an integer index rather than current density for an eigenstate)

x_{j} = \frac{x_{o}}{y} + (2 j - 3) \frac{L}{2},

(64)

L = x_{2} - x_{1} = \frac{2 x_{o}}{y} \sqrt{1 - y^{2}},

(65)

κ = \frac{2}{ℏ} {[\frac{m Q}{x_{o}} \frac{(1 - y)}{y}]}^{1 / 2},

(66)

where the energy dependence of

x_{j}, L

⁠, and

κ

is suppressed but implicit in the definition of

y

on which they depend. The area factor for the SN barrier can be simplified using

x \equiv [(x_{2} + x_{1})

+ (x_{2} - x_{1}) \cos φ] / 2

in Eq. to reveal the explicit

y

-dependence

\begin{aligned} σ_{s n} (y) & = \int_{x_{1}}^{x_{2}} {[\frac{x_{o} (x_{2} - x) (x - x_{1})}{x (x_{2} - x_{o}) (x_{o} - x_{1})}]}^{1 / 2} d x \\ = \frac{\sqrt{1 + y}}{2 \sqrt{2}} \int_{0}^{π} \frac{\sin^{2} φ d φ}{{[1 + \sqrt{1 - y^{2}} \cos φ]}^{1 / 2}}, \end{aligned}

(67)

for which analytic integration shows that

σ_{s n} (0) = σ_{2} = 2 / 3

(the triangular barrier) and

σ_{s n} (1) = σ_{\cap} = π / 4

(the parabolic barrier), as shown in Fig. 5. Combining Eqs. , , and into

θ = 2 σ κ L

gives

\begin{matrix} θ_{s n} (y) = \frac{2 \sqrt{2 m Q x_{o}}}{ℏ} \frac{(1 - y^{2})}{y^{3 / 2}} F [p (y)], \\ F (p) \equiv \int_{0}^{π} \frac{\sin^{2} φ}{\sqrt{1 + p \cos φ}} d φ, \\ p (y) \equiv \sqrt{1 - y^{2}}, \end{matrix}

(68)

where

x_{o} = \sqrt{Q / F}

is used in preference to

F

⁠. The value of

F [p (y)]

at

p = (0, 1)

[corresponding to

y = (1, 0)

] is analytic and given by

F (1) = 4 \sqrt{2} / 3

⁠, and

F (0) = π / 4

⁠. Similar to the behavior of

σ_{s n} (y)

⁠, direct evaluation shows that a weak convexity for

F [p (y)]

vs

(1 - y) / (1 + y)

exists, and so an accurate approximation is⁹²

\begin{matrix} F_{n} [p (y)] \equiv \sum_{j = 0}^{n} C_{j} {[\frac{1 - y}{1 + y}]}^{j}, \\ F_{n}^{'} [p (y)] \equiv \sum_{j = 1}^{n} j C_{j} {[\frac{1 - y}{1 + y}]}^{j - 1}, \end{matrix}

(69)

where representative cases are shown in Table III, and where

F_{n}^{'} [p (y)]

is defined for later usage, corresponding to the derivative of

F_{n}

with respect to the polynomial argument

[(1 - y) / (1 + y)]

⁠. Even for the lowest

n = 2

case, as Fig. 6 shows, the absolute accuracy is better than

0.11 %

so that the difference between

F

and

F_{n}

is not easily discernible: consequently, for standard numerical evaluations,

F_{2} [p (y)]

is sufficiently accurate.

FIG. 5.

View large Download slide

Comparison of the Schottky–Nordheim barrier $σ (y)$ exactly evaluated by Eq. (67) and by a quadratic approximation in $x (y)$ based on three points calculated at $x = [0, 1 / 2, 1]$ for $x (y) = (1 - y) / (1 + y)$ ⁠.

FIG. 6.

A comparison of F [ p ( y ) ] − F ( 0 ), where p ( y ) = 1 − y 2 and F ( 0 ) = 4 2 / 3, to its approximation by Eq. (69) for n = ( 2 , 4 , 6 ): (red, green, and blue dashed lines overlapping ∙ points). To show the accuracy, error plots defined by E r r n = 100 % × [ 1 − ( F n / F ) ] are shown, scaled by a factor of 3 × to show detail: for n = 2, | E r r 2 | < 0.11 %.

View large Download slide

A comparison of $F [p (y)] - F (0)$ ⁠, where $p (y) = \sqrt{1 - y^{2}}$ and $F (0) = 4 \sqrt{2} / 3$ ⁠, to its approximation by Eq. (69) for $n = (2, 4, 6)$ ⁠: (red, green, and blue dashed lines overlapping $∙$ points). To show the accuracy, error plots defined by $E r r_{n} = 100 % \times [1 - (F_{n} / F)]$ are shown, scaled by a factor of $3 \times$ to show detail: for $n = 2$ ⁠, $| E r r_{2} | < 0.11 %$ ⁠.

TABLE III.

Coefficients for the polynomial approximation to $F_{n} [p (y)]$ in Eq. (69). In all cases, the exact values C₀ = π/4 and $\sum_{j = 1}^{n} C_{j} = 4 \sqrt{2} / 3$ can be used.

C_j	n = 2	n = 4	n = 6
C₀	1.570 80	1.570 80	1.570 80
C₁	0.580 164	0.591 405	0.590 150
C₂	−0.265 342	−0.335 370	−0.334 254
C₃	…	0.131 396	0.178 278
C₄	…	−0.072 608 7	−0.250 222
C₅	…	…	0.226 839
C₆	…	…	−0.095 968 7
$\sum C_{j}$	1.885 62	1.885 62	1.885 62

C_j	n = 2	n = 4	n = 6
C₀	1.570 80	1.570 80	1.570 80
C₁	0.580 164	0.591 405	0.590 150
C₂	−0.265 342	−0.335 370	−0.334 254
C₃	…	0.131 396	0.178 278
C₄	…	−0.072 608 7	−0.250 222
C₅	…	…	0.226 839
C₆	…	…	−0.095 968 7
$\sum C_{j}$	1.885 62	1.885 62	1.885 62

Although applied to the SN barrier, the method of expressing $θ (y)$ in terms of $y$ -dependent coefficients from $2 σ κ L$ and a generally smoothly varying polynomial originating from $σ (y)$ is, in fact, a general method: variations on it can be used for more complex potentials characteristic of space charge, interface barriers, and curved trajectories from hyperbolic emitters.^86,88

6. Modified Schottky–Nordheim barriers

The SN barrier is based on a charge outside the surface of a planar metallic conductor. The form of the image charge can be trivially modified by

Q \to (K_{s} - 1) Q / (K_{s} + 1)

⁠, where

K_{s}

is the dielectric constant of the semiconductor.⁹³ However, the more consequential modifications to the SN are a consequence of

quantum deviations from the classical image charge,^94–96
changes to the barrier itself due to nanometric curvature and space charge effects that can be approximated by so-called “quadratic” or “parabolic” extensions^86,97,98 to model conical and ellipsoidal emitters, and
modifications to the form of the image charge for a curved surface.^99,100

The quadratically modified barriers can be dealt with computationally by adding a term

γ x^{3}

to

V (x)

⁠, resulting in a requirement for finding the roots of a cubic equation so as to rapidly evaluate the Gamow factor using shape factor methods. The most difficult case, however, is accommodating the modifications to the image charge term. The modifications are well-known to be

\frac{Q}{x} \to \frac{4 a Q}{r^{2} - a^{2}} \approx \frac{Q}{x} - \frac{Q}{2 a} + \frac{Q x}{4 a^{2}},

(70)

for a charge outside a sphere of radius

a

⁠, where

r = x + a

⁠. It is seen that defining effective parameters

Φ_{e f f} = Φ + Q / 2 a

and

F_{e f f} = F - Q / 4 a^{2}

allow the retention of the SN formalism both in the canonical emission equations and the shape factor method, and so the canonical equations or their recent extensions can be applied to small patches tiling the surface to evaluate total emitted current from a protrusion.¹⁹ A similar comment applies to changes due to exchange-correlation effects that offset the origin by

x_{o}

and make the planar image charge resemble

Q / (x + x_{o})

⁠.

More difficult, however, is the case of an image charge outside oft-used ellipsoidal or hyperbolic shape, so-called prolate spheroidal models, because there is no single radius of curvature as for the sphere. The definition of an effective radius is more nuanced. In that case, a “mean” curvature can be found: if the surface height is defined by a function

z_{s} = g (x, y)

⁠, and

g_{s} = \partial_{s} g, g_{s s^{'}} = \partial_{s} \partial_{s^{'}} g

for

(s, s^{'}) \in (x, y)

⁠, then Kyritsakis¹⁰⁰ has shown that the effective radius of curvature to use is

\frac{1}{a_{e f f}} = - \frac{(1 + g_{y}^{2}) g_{x x} - 2 g_{x} g_{y} g_{x y} + (1 + g_{x}^{2}) g_{y y}}{2 {(1 + g_{x}^{2} + g_{y}^{2})}^{3 / 2}},

(71)

and is used to specify to the aforementioned parabolic factor

γ

for curvature-affected barriers.¹⁰⁰ For prolate spheroidal coordinates,^8,101,102 the surface can be analytically defined and the gradients evaluated.

For modeling purposes, a curved surface image charge and a prolate spheroidal barrier both introduce a positive quadratic term $γ x^{2}$ ⁠, which causes the width of the barrier $L (E)$ to enlarge for a given energy $E$ and the shape factor $σ (y)$ to reduce. That will cause the current density to decrease for the canonical equations. For field emission in particular, the magnitude of the surface field declines down the sides of the emitter, which causes the tunneling barrier width at $E$ [the factor $L (E)$ in $θ (E)$ ] to enlarge and thereby causes $D (E)$ to become smaller.⁸⁶ The effect contributes to a focusing behavior termed “self-focusing” by Kyritsakis et al.⁹⁸

III. CURRENT DENSITY FORMULATIONS

The current density

J (F, T)

for general barriers can be numerically evaluated, but analytical forms of great utility are possible by introducing simple approximations.^8,103 Let the barrier height be

V_{o} = μ + ϕ

⁠, where

μ

is the Fermi level, and although

ϕ

is motivated by the SN barrier, the methods to follow are general. First, rewrite Eq. in terms of an integration over energy

E

rather than wave number

k

using the assumption that energy is parabolic in momentum,

E = ℏ^{2} k^{2} / 2 m

⁠, and no electrons are incident from the right,

J (F, T) = \frac{q}{2 π ℏ} \int_{0}^{\infty} D (E) f (E) d E .

(72)

Second, linearize the Gamow factor about a suitable expansion point

E_{m}

such that

\begin{matrix} θ [y (E)] \approx θ (y_{m}) - β_{F} (y_{m}) (E_{m} - E), \\ β_{F} (y) = - \frac{d}{d E} θ [y (E)] = - \frac{y^{2}}{\sqrt{4 Q F}} \frac{d}{d y} θ (y), \end{matrix}

(73)

where

β_{F}

will be referred to as the energy slope factor for field and

y_{m} \equiv y (E_{m})

⁠. Equation shall be called the “Linear-

θ (E)

” (⁠

L θ

⁠) approximation below. The choice of

E_{m}

depends on where the maximum of

D (E) f (E)

in the integrand is located, but two asymptotic limits in

F

result in the thermal and field approximations leading to the RLD (thermal) and FN (field) equations and so consider those limits exclusively here.

A. Thermal dominated

1. Temperature determination

The temperature

T

in the supply function is that of the electron gas, assumed to be the same as that of the lattice of atoms making up the cathode material. Often, however, the temperature of the cathode is determined by some standard method such as a pyrometer, where gray body emission is compared to black body emission to infer

T

from the (pyrometer) brightness temperature

T_{b}

by means of equating gray body radiation spectrum to black body (Planck’s law)^1,104 and eliminating common factors to obtain

\frac{t r ϵ}{\exp (ℏ ω / k_{B} T) - 1} = \frac{1}{\exp (ℏ ω / k_{B} T_{b}) - 1},

(74)

where

ω = 2 π c / λ

⁠, and

(t, r, ϵ)

are the transmissivity of the viewport, reflectivity of the mirror, and the emissivity of the cathode, respectively. Pyrometers typically chose one or more wavelengths from which to make a measurement (e.g., dual wavelength pyrometers might choose 700 and 800 nm). The solution of Eq. for

T (λ, T_{b})

is shown in Fig. 7.

FIG. 7.

View large Download slide

Solution of Eq. (74) for actual temperature $T$ as a function of wavelength $λ$ and brightness temperature $T_{b}$ (shown in legend) for typical dispenser cathode values¹⁰⁴ of $(t, r, ϵ) = (0.93, 0.76, 0.52)$ ⁠, respectively.

2. Richardson–Laue–Dushman equation

For thermal emission,

β_{F} (y)

is large even near the barrier maximum where

θ (1) = 0

identically, so that the Kemble approximation means that

D (E)

resembles a Heaviside step function in the region where

f (E)

is non-negligible, making

E_{m} \approx μ + ϕ

⁠. It is, therefore, sufficient to keep the approximation of a step function, giving for the thermal limit

\begin{aligned} J_{T} (T) & \approx \frac{q}{2 π ℏ} \int_{μ + ϕ}^{\infty} \frac{m}{π β_{T} ℏ^{2}} e^{β_{T} (μ + ϕ - E)} d E = \frac{q m}{2 π^{2} β_{T}^{2} ℏ^{3}} e^{- β_{T} ϕ}, \end{aligned}

(75)

where

β_{T} = 1 / k_{B} T

⁠. Introducing

A_{R L D} = \frac{q m k_{B}^{2}}{2 π^{2} ℏ^{3}},

(76)

numerically equal to

7.500 64 \times 10^{- 9}

(q/fs nm

^{2}

K

^{2}

⁠) in (efnq) or 120.1735 (A/cm

^{2}

K

^{2}

⁠) in SI units results in the conventional RLD equation

J_{R L D} = A_{R L D} T^{2} \exp [- ϕ / k_{B} T],

(77)

where the field dependence of

ϕ = Φ - \sqrt{4 Q F}

is implicit. For

ϕ = 2

eV and

T = 1273

K, then

J_{R L D} = 2.352 45

A/cm

^{2}

⁠. When

D (E)

cannot be treated as a step function, a common recourse is to append a factor

(1 - r)

to

J_{R L D}

and interpret

r

as a mean quantum mechanical reflection term,¹⁰⁵ e.g.,

r \approx 0.05

or so, but under such conditions, recourse to the General Thermal-Field (GTF) limit in Sec. III C 1 is preferable.

It is evident that the simplicity of Eq. (77) derives in part from the transverse directions in three dimensions (3D) being infinite in extent, thereby allowing the integrations in Eq. (42) to be performed. However, thermionic emission is also possible from two-dimensional (2D) materials like graphene. Treating electrons in graphene like a Dirac fermion causes the number of states between an energy $E_{p}$ and $E_{p} + d E_{p}$ to be given by $2 E_{p} d E_{p} / π {(ℏ v_{f})}^{2}$ ⁠, where $v_{f}$ is the massless Dirac fermion velocity in the graphene. As a consequence, whereas the RLD equation entails $\ln [J_{R L D} / T^{2}]$ is linear in $1 / T$ ⁠, for graphene, $\ln [J_{2 D} / T^{3}]$ is linear in $1 / T$ for thermal emission normal to the sheet.¹⁰⁶

B. Field dominated

1. Simple barriers

For tunneling emission in the low temperature limit, Eq. results in the approximation

\begin{aligned} J_{F} (F) & = \frac{q m e^{- θ (y_{m})}}{2 π^{2} ℏ^{3} β_{F} {(y_{m})}^{2}} \int_{0}^{β_{F} μ} s e^{- s} d s \approx \frac{q m e^{- θ (y_{m})}}{2 π^{2} ℏ^{3} β_{F} {(y_{m})}^{2}} \end{aligned},

(78)

in the limit that

β_{F} (y_{m}) μ ≫ 1

⁠. Observe that now both

θ

and

β_{F}

use

y

as their arguments rather than energy and that

y_{m}

is defined by

E_{m}

⁠: these observations depart from the conventions of prior work [where energy was the argument of

θ

and

β_{F}

⁠, and

y_{o} \equiv y (μ)

referred to its value for the SN barrier at

E_{m} \approx μ

]. The change is because the present methods are developed to apply to general barriers as well.

The two simplest and common cases for tunneling emission are the triangular barrier of the original Fowler–Nordheim equation and the SN barrier forms that led to its refinement. The rectangular barrier case is trivial and need not be explicitly considered.

For the triangular barrier,

σ_{▹} = 2 / 3

⁠, giving

\begin{matrix} θ_{▹} (E) = \frac{4}{3 ℏ F} \sqrt{2 m {(μ + Φ - E)}^{3}}, \\ β_{F} (E) = \frac{2}{ℏ F} \sqrt{2 m (μ + Φ - E)}, \end{matrix}

(79)

where the forms are in terms of energy because

V_{o} = μ + Φ

and so

V_{o} = V (x_{o})

because

Q = 0

in Eq. for the triangular barrier. The analytic forms become for

E_{m} \approx μ

⁠,

\begin{matrix} θ_{▹} (μ) = \frac{4}{3 ℏ F} \sqrt{2 m Φ^{3}}, \\ β_{F} (μ) = \frac{2}{ℏ F} \sqrt{2 m Φ} . \end{matrix}

(80)

For a work function of

Φ = 4.5

eV and

μ = 7

eV, then

y (μ) = 0.754 238

⁠,

θ_{▹} (μ) = 8.150 91

⁠, and

β_{F} (μ) = 2.716 97

eV

^{- 1}

for

F = 8

eV/nm so that

β_{F} (μ) μ = 1 9.018 79

⁠. As a consequence,

e^{- β_{F} μ} ≪ 1

⁠, and so Eq. is well approximated by

J_{F} \approx \frac{q}{16 π^{2} ℏ Φ} F^{2} \exp [- \frac{4 \sqrt{2 m Φ^{3}}}{3 ℏ F}] .

(81)

It follows that

\ln [J_{F} / F^{2}]

is linear in

1 / F

⁠, analogous to the behavior of

\ln [J_{T} / T^{2}]

being linear in

1 / T

from Eq. . As with thermionic emission from 2D materials, changes in the energy dispersion relations for two-dimensional materials like graphene considered before, so too will those changes affect the field emission equation by rendering

\ln [J_{2 D} / F^{3}]

linear in

1 / F

for field emission.¹⁰⁷

The right hand side of Eq. is almost the equation that Fowler and Nordheim⁸⁷ derived: because they effectively used asymptotic relations for Airy functions to solve Schrödinger’s equation for the triangular barrier, it can be shown [Ref. 8, Eq. (13.9)] that the triangular barrier FN equation is a minor modification given by

J_{F N} (F) = \frac{4 \sqrt{μ Φ}}{μ + Φ} J_{F} .

(82)

The extension of

J_{F}

to the SN barrier to include Schottky factor effects causes

σ

to no longer be constant as a function of

y (E)

⁠. As a result, although Eq. corresponds to the triangular barrier model of Fowler and Nordheim,⁸⁷ it differs from their later model based on the SN barrier, for which Murphy and Good (MG)²⁴ give the most commonly used form. The MG form will be considered after completing the analysis of the SN barrier.

Last, at high fields, the SN barrier above the Fermi level takes on a parabolic shape: in the limit that it is parabolic, then the parabolic approximations are

\begin{matrix} θ_{\cap} (μ + ϕ) = 0, \\ β_{F} (μ + ϕ) = \frac{π}{ℏ F} \sqrt{m Φ y (μ)}, \end{matrix}

(83)

where

y (μ) = \sqrt{4 Q F} / Φ

and

ϕ = [1 - y (μ)] Φ

describe conditions near the barrier maximum.⁵

Some points concerning $J_{T}$ and $J_{F}$ for the simple barriers require emphasis. First, for Eq. (75), $J_{T}$ appears to be only a function of temperature $T$ ⁠, but if $ϕ$ has a field dependence, then it is implicitly a function field as well. Second, because the zero temperature limit of the supply function was used for Eq. (78), $J_{F}$ only has a field dependence, but under warm conditions, the supply function introduces a temperature dependence. When both are generalized in the GTF limit, both become functions of $T$ and $F$ directly.

2. Schottky–Nordheim barrier

The derivative with respect to energy

E

of

θ [y (E)]

can be taken for the SN barrier using Eq. , which allows derivatives in terms of

E

to be re-expressed in terms of

y (E)

⁠. The derivative has two parts: the first finds the derivative of

\sqrt{1 - y^{2}} / y^{3 / 2}

and the second of

F [p (y)]

⁠. They are

\begin{aligned} y^{2} \frac{d}{d y} [\frac{\sqrt{1 - y^{2}}}{y^{3 / 2}}] & = - \frac{3 - y^{2}}{2 \sqrt{y (1 - y^{2})}}, \\ y^{2} \frac{d}{d y} F_{n} [p (y)] & = - \frac{2 y^{2}}{{(1 + y)}^{2}} F_{n}^{'} [p (y)], \end{aligned}

where

F^{'}

is defined in Eq. . Combining these results gives

β_{F} (y) = \frac{1}{ℏ} \sqrt{\frac{m x_{o}^{3}}{2 Q y}} [(3 + y^{2}) F_{n} + 4 y (\frac{1 - y}{1 + y}) F_{n}^{'}],

(84)

where for visual clarity, the arguments of

F_{n} [p (y)]

and

F_{n}^{'} [p (y)]

are suppressed. The factor in brackets is equal to

2 π

when

y = 1

⁠. From the forms of

θ_{s n} (y)

[Eq. ] and

β_{F} (y)

[Eq. ], it is clear that both are the product of a coefficient that varies only with

F

multiplying a factor depending only on

y

⁠. Consequently, the formulation can be rapidly implemented numerically, as the values of

C_{j}

are unchanging.

The leading order approximation to

J_{F} (F)

for the Schottky–Nordheim barrier and

y_{m} = y (E_{m})

is then

J_{F} (F) = \frac{q m}{2 π^{2} ℏ^{3} β_{F} {(y_{m})}^{2}} \exp [- θ (y_{m})],

(85)

where

E_{m}

must be determined from the maximum of

D (E) f (E)

in the current density integrand, which in general is found by numerical means. The convexity of

θ [y (E)]

as a function of

[(1 - y) / (1 + y)]

and, therefore, of

1 / y (E) = [μ + Φ - E] / \sqrt{4 Q F}

in Figs. 8 and 9 reveal that where the expansion point

E_{m}

is taken in Eq. can be consequential if the temperature is such that the peak of

D (E) f (E)

moves into the region between

μ

and

μ + ϕ

⁠, as occurs during the transition from field to thermal emission (Fig. 9).

FIG. 8.

View large Download slide

The Gamow factor $θ_{s n} [y (E)]$ [Eq. (68)] at various $F$ and its ratio to $θ [y (μ)]$ [Eq. (10)]. Dashed lines correspond to FN approximation of Eq. (87). Yellow dot is unity. $μ = 7$ eV, $Φ = 4.5$ eV.

FIG. 9.

View large Download slide

The energy slope factor $β_{F} (y)$ [Eq. (84)] and its ratio to $β_{F} (1)$ at shown fields. Gray dots show $C_{F N} / β_{F} (1)$ using Eqs. (87) for the same fields. $μ = 7$ eV, $Φ = 4.5$ eV.

3. Fowler–Nordheim equation

In the case of zero temperature, though,

f (E)

vanishes for

E > μ

⁠, and so

E_{m}

is always below the Fermi level. Fowler and Nordheim made the convenient assumption that

E_{m} = μ

is at the Fermi level, and that the Fermi level was well below the barrier maximum

μ + ϕ

⁠. Their choice of

y_{m} \to y_{o} = y (μ) = \sqrt{4 Q F} / Φ

= 2 Q / x_{o} Φ

in

J_{F} (F)

results in the most common form used of the Fowler–Nordheim equation as rendered by Murphy and Good,^24,25 although a modern updating of the theory by Forbes and Deane¹⁰⁸ provides the most elegant representation of the elliptical integral function terms

v (y)

and

t (y)

that Murphy and Good exploit. The Forbes–Deane (FD) approximations are given by

\begin{matrix} v (y) \approx 1 - \frac{y^{2}}{3} [3 - \ln (y)], \\ t (y) \approx 1 + \frac{y^{2}}{9} [1 - \ln (y)], \end{matrix}

(86)

and are referred to as “Schottky–Nordheim barrier functions” by them. The linearization of the Gamow factor about the expansion point

E = μ

results in

\begin{matrix} θ_{f n} (E) = B_{f n} + C_{f n} (μ - E), \\ B_{f n} \equiv \frac{4}{3 ℏ F} \sqrt{2 m Φ^{3}} v (y_{o}), \\ C_{f n} \equiv \frac{2}{ℏ F} \sqrt{2 m Φ} t (y_{o}), \end{matrix}

(87)

where

y_{o} = y (μ)

⁠,

B_{F N} \leftrightarrow θ (y_{o})

⁠, and

C_{F N} \leftrightarrow β_{F} [y (μ)]

when compared to Eq. and used in Eq. . For typical metal-like values (⁠

μ = 7

eV,

Φ = 4.5

eV,

F = 8

eV/nm), then for

y_{o} = 0.754 238

⁠, Eq. gives

v (y_{o}) = 0.377 642

and

t (y_{o}) = 1.081 04

⁠.

From the definitions of

θ [y (μ)]

and

β_{F} [y (μ)]

⁠, it follows directly that the shape factor representation of the elliptical integral functions is given by

\begin{matrix} v (y) = \frac{3 \sqrt{2}}{8} (1 - y^{2}) F_{n} [p (y)], \\ t (y) = \frac{\sqrt{2}}{8} [(3 + y^{2}) F_{n} + 4 y (\frac{1 - y}{1 + y}) F_{n}^{'}], \end{matrix}

(88)

where the arguments of

F_{n} [p (y)]

and

F_{n}^{'} [p (y)]

are suppressed for visual clarity. The accuracy of Eq. , even for

n = 2

⁠, is significantly superior for the range of

y

typical of field emission conditions

(y ≳ 0.5)

to Eq. as shown in Fig. 10, but an error of 0.33% is generally difficult to distinguish visually. As far as numerical convenience is the relevant metric, they are comparable with FD having an advantage, but there are other considerations. First, the shape factor approach is generalizable to other potential barriers⁸⁶ but the usage of the Schottky–Nordheim barrier functions is restricted to a triangular barrier modified by an image charge potential. Second, when used in the MG formulation of the FN equation, the FD approximations yield a convenient form given by¹⁰⁹^,

\begin{matrix} J_{M G} (F, T = 0) \equiv A_{M G} {(\frac{F}{B_{o}})}^{2 - ν} \exp (- \frac{B_{o}}{F}), \\ A_{M G} \equiv \frac{A_{o} B_{o}^{2}}{Φ t_{o}^{2}} {(\frac{Φ^{2} e^{6}}{4 Q B_{o}})}^{ν}, \end{matrix}

(89)

where

A_{o} \equiv q / 16 π^{2} ℏ = 0.009 620 87

(q/eV fs),

ν = (8 Q / 9 ℏ)

\sqrt{2 m / Φ}

⁠,

B_{o} \equiv (4 / 3 ℏ) \sqrt{2 m Φ^{3}}

⁠, and

t_{o} \equiv t (y_{o})

⁠. Observe that

1 < t (y) < 10 / 9

exhibits a weak variation with field. Representative values for

Φ = 4.5

eV are

B_{o} = 65.2073

eV/nm,

ν = 0.772 808

⁠, and

(Φ^{2} e^{6} / 4 Q B_{o}) = 87.0050

⁠. To reflect its genesis, Eq. should be described as the FN equation with image charge methods of Murphy and Good using the approximations of Forbes and Deane, but it is occasionally and not entirely correctly referred to as simply the Fowler–Nordheim equation. In deference to its legacy, Eq. shall be denoted by the MG subscript, based as it is on a simplification of the Murphy and Good equation. Its outstanding convenience is two-fold. It is convenient because, first,

{(Φ^{2} e^{6} / 4 Q B_{o})}^{ν}

describes how much the image charge modification increases the current in excess of the triangular barrier (or

31.5424

for

Φ = 4.5

eV), and second, because it allows factors that affect total current (current density integrated over emission area) such as statistics^110,111 and non-linear barriers⁸⁶ to simply appear as modifications to

{(F / B_{o})}^{2 - ν}

in the current density equation itself (apart from other multiplying coefficients such as notional emission area^8,111).

FIG. 10.

View large Download slide

The accuracy of the equations for $v (y)$ and $t (y)$ using the Forbes–Deane approximation [Eq. (86)] compared to the shape factor approach [Eq. (88)] for $n = (2, 4)$ ⁠. For field emission from metals, $y$ ≳ $0.5$ ⁠.

The extension of

J_{M G} (F, T)

to include a temperature dependence was achieved by Murphy and Good²⁴ by appending a term that arose when

β_{T}

was finite. Although it shall be replaced by the GTF formulation below, for completeness (see Sec. 24.3.2 of Ref. 8 or Ref. 25), it results in

J_{M G} (F, T) = J_{M G} (F, 0) \frac{π k_{B} T C_{f n}}{\sin (π k_{B} T C_{f n})},

(90)

where Eq. for

C_{f n}

is required. The correction is singular when

C_{f n} = 1 / k_{B} T

and so the temperature-dependent addition cannot track the entry of current density into the thermal-field regime from the field regime side, and so it is not considered further. Below,

J_{M G} (F)

shall refer to

J_{M G} (F, 0)

exclusively.

A comparison of $J_{M G}$ to $J_{R L D}$ reveals the domains where each is appropriate: where they are equal divides those regions and constitutes the basis of nexus theory.^112–114 The comparison as a function of field, with $J_{R L D}$ evaluated for various temperatures, is shown in Fig. 11. Clearly, there are regimes where field emission dominates thermal emission. The abrupt transition from thermal emission behavior to field emission shown in that figure accounts for similar features seen in $I (V)$ data from thermionic emitters by Geittner et al.¹¹⁵ (see also Chap. 8 in Ref. 116), where field emission contributions from protrusions on a surface contribute to the thermal emission from a Ba-dispenser cathode. A second example using the data of Dyke et al.³⁴ to show the important impact of temperature on field emission, must await the development of the $r$ GTF method in Sec. III C 2.

FIG. 11.

View large Download slide

Comparison of $J_{M G} (F, 0)$ [Eq. (89), black dashed] to various temperatures of $J_{R L D} (F, T)$ [Eq. (77)] for the metal-like parameters $Φ = 4.5$ eV. A standard $T$ vs $J_{R L D} (F, T)$ plot would use the data points along a vertical line for a given $F$ ⁠.

C. Thermal-field emission

A rapid appreciation of the complexity thermal-field emission introduces is made possible using the triangular barrier transmission probability of Eq. (79), as it is entirely analytic. An extension to the SN barrier is straightforward although its numerical implementation is more nuanced. In either case, however, the finding of the maximum location $E_{m}$ of the product $D (E) f (E)$ from which $y (E_{m})$ is determined proceeds analogously. Consider two complementary cases: first, hold the field fixed at a representative value and increase the applied temperature as in Fig. 12, and second, hold the temperature fixed and increase the applied field as in Fig. 13. In both cases, parameters are set so that $μ + Φ \geq E_{m} \geq μ$ ⁠, that is, between the thermal and field regimes. Although on a log-plot, it is evident that the broadness of $D (E) f (E)$ is substantially larger than for either thermal or field emission. Note that the small perturbations near $E \approx μ + Φ \equiv V_{o}$ in $D (E)$ are due to both $θ (V_{o})$ and $β_{F} (V_{o})$ vanishing, as per Eq. (79), which causes a vaguely discernable double hump to appear in $D (E) f (E)$ there.

FIG. 12.

View large Download slide

Gray line: $D (E)$ for $μ = 7$ eV, $Φ = 2$ eV, $F = 0.5$ eV/nm. Thin color lines: $f (E)$ ⁠; Thick color lines: $D (E) f (E)$ ⁠, both labeled by $T$ (K). The black bullet (⁠ $∙$ ⁠) is $D (E_{m}) f (E_{m})$ ⁠. The blue region is for $E < μ$ ⁠; the orange region is for $E > μ + Φ$ ⁠.

FIG. 13.

View large Download slide

Gray line: $f (E)$ for $μ = 7$ eV, $Φ = 4$ eV, $T = 1500$ K. Thin color lines: triangular barrier $D (E)$ ⁠; Thick lines: $D (E) f (E)$ ⁠, both labeled by $F$ (eV/nm). The black bullet (⁠ $∙$ ⁠): $D (E_{m}) f (E_{m})$ ⁠. Ther blue region is for $E < μ$ ⁠; the orange region is for $E > μ + Φ$ ⁠.

The behavior of

D (E) f (E)

in the thermal-field gap shown in Fig. 14 is qualitatively understood as follows. In that region,

f (E) \approx \exp [- β_{T} (E - μ)]

⁠, whereas

D (E) \approx \exp [- θ (E_{m}) - β_{F} (E_{m})

(E_{m} - E)]

⁠. The maximum occurs where

\partial_{E} \ln [D (E) f (E)] = 0

⁠, which, therefore, entails that

β_{T} - β_{F} (E_{m}) = 0

⁠. The simplicity of the triangular barrier allows this equation to be solved for

E_{m}

as a function of field using Eq. , giving the condition for thermal-field emission for the triangular barrier,

E_{m} \approx μ + Φ - \frac{1}{2 m} {[\frac{ℏ F}{2 k_{B} T}]}^{2} .

(91)

The behavior of

n (F, T) \equiv \frac{β_{T}}{β_{F} (E_{m})},

(92)

is, therefore, seen to partition the thermal (⁠

n ≲ 0.95

⁠), field (⁠

n ≳ 1.05

⁠), and thermal-field (⁠

0.95 < n < 1.05

⁠) regimes. Introducing

z = β_{F} (E - μ)

and using

n (F, T)

in favor of the energy slope factors results in the compact thermal-field current density. Although the analysis was for a triangular barrier, the same methodology holds for more complex barriers like the SN barrier: all that is required is that a means to calculate both

θ (E_{m})

and

β_{F} (E_{m})

exists, analytically if possible, numerically if not. The shape factor approach to evaluate

θ (E) = 2 σ [y (E)] κ (E) L (E)

provides a fast method to determine the Gamow factor and its derivatives.⁸⁶

FIG. 14.

View large Download slide

Contour plot of current density integrand $D (E) f (E)$ for $Φ = 4$ eV and $T = 1500$ K. Open bullets (⁠ $°$ ⁠) are numerically determined location of $E_{m}$ ⁠: the data points are represented again in Fig. 15. Thick yellow line is Eq. (91) for $F$ in the thermal-field regime.

1. Theory

Return now to Eq. and write it in a form that incorporates Eq. (supply function), Eq. (transmission probability), and Eq. (dimensional coefficient) but now using the parameters appropriate for thermal and field emission. It becomes

J = A_{R L D} T^{2} \int_{0}^{\infty} \frac{\ln {1 + \exp [β_{T} (μ - E)]}}{1 + \exp [θ (E)]} β_{T} d E,

(93)

where the integral has been put into a form that is purposely dimensionless. It is desired to recast the integral so that

θ

becomes the integration term: if, to a good approximation,

θ (E) \approx θ (E_{m}) + β_{F} (E_{m}) (E - E_{m})

⁠, then

\begin{matrix} θ (E) \approx β_{F} (E_{m}) (E_{o} - E), \\ E_{o} \equiv \frac{θ (E_{m})}{β_{F} (E_{m})} + E_{m} . \end{matrix}

(94)

Now, letting

z \equiv β_{F} (E_{m}) (E_{o} - E)

results in

\begin{matrix} β_{T} (μ - E) = n β_{F} [- (E_{o} - μ) + (E_{o} - E)] \equiv n (- s + z) . \end{matrix}

(95)

These replacements and the approximation

β_{F} E_{o} ≫ 1

allow Eq. to be rewritten in terms of the new parameters

n (F, T)

= β_{T} / β_{F} (E_{m})

⁠,

s (F, T) = β_{F} (E_{m}) (E_{o} - μ)

⁠, and

E_{o} = E_{m}

+ θ (E_{m}) / β_{F} (E_{m})

as

\begin{matrix} J (F, T) = A_{R L D} T^{2} N (n, s), \\ N (n, s) \equiv n \int_{- \infty}^{\infty} \frac{\ln [1 + e^{n (z - s)}]}{1 + e^{z}} d z . \end{matrix}

(96)

Although

N (n, s)

can be found through numerical integration, it is an unpleasant combination of factors that separately can be large, e.g., if

E_{m} \approx μ

⁠, then

s \approx θ (μ)

⁠. For a triangular barrier for

F = 6

eV/nm and

Φ = 4.5

eV,

θ_{▹} (μ) \approx 4 \sqrt{2 m Φ^{3}} / 3 ℏ F = 10.8679

⁠; a Schottky barrier reduces that by about a factor of 2. As

n

can range from small to large, care is then required to monitor the integration, undercutting the utility of Eq. in simulations. A laborious expansion of the functions in the integrand followed by a term-by-term integration and summation of related terms produces a hierarchy of contributions^8,103 for which the dominant ones can be grouped, allowing for an analytical approximation of reasonable accuracy with rapid execution speed given by

\begin{matrix} N (n, s) \approx e^{- s} n^{2} Σ (\frac{1}{n}) + e^{- n s} Σ (n), \\ Σ (x) \approx \frac{1 + x^{2}}{1 - x^{2}} - a x^{2} - b x^{4}, \end{matrix}

(97)

where

a = ζ (2) - 2 = - 0.355 066

and

b = (7 ζ (4) - 8) / 4

= - 0.105 934

⁠, where

ζ (x)

is Riemann’s Zeta function.¹¹⁷ A measure of the quality of the approximation is suggested by observing that in the limit that

n \to 1

⁠, the approximation approaches

N (1, s) \to (s + 1) e^{- s}

⁠, which can be shown to be the exact value.¹¹⁸ For the SN barrier, the relationship of

J (F, T)

from Eq. to the prior canonical equations is

lim_{n \to \infty} A_{R L D} T^{2} N (n, s) \approx J_{F N} (F),

(98)

lim_{n \to 0} A_{R L D} T^{2} N (n, s) \approx J_{R L D} (F) .

(99)

Although close, the equality is not exact because of the contributions of

Σ (x)

and the departure of

E_{m}

from

μ

for field emission or

D (μ + ϕ) < 1

for thermal emission. It is seen that the

n^{2} e^{- s}

part of

N (n, s)

is the “field-like” contribution, and the

e^{- n s}

part is the “thermal-like” contribution. The “thermal-field” regime is for

n \approx 1

⁠, where the boundaries of

0.95 \leq n \leq 1.05

are useful as in Fig. 15. Further,

Σ (x)

has a singularity at

x \to 1

⁠, and so both

Σ (n)

and

Σ (1 / n)

must be retained as

n

approaches unity (i.e., both the thermal and field contributions must be combined) to insure the singularity is mitigated. Finally,

Σ (1 / n)

gives the thermal contributions to field emission, and

Σ (n)

gives the field contributions to thermal emission. A comparison of

J (F, T)

from Eq. is shown in Fig. 16: the broad width of the plateau of the integrand in the thermal-field region for

n \approx 1

is a consequence of the approximation that

β_{F} (E_{m}) (E - μ) \leftrightarrow n (z - s)

is made to be linear in

E

⁠. As is apparent from Fig. 8, however,

θ (E)

is only approximately linear in

E

⁠, and because

D (E)

depends on

e^{θ}

⁠, the non-linearities can have an outsized effect.

FIG. 15.

View large Download slide

$F$ vs $n (F, T)$ for $Φ = 4$ eV and $T = 1500$ K. Open bullets (⁠ $°$ ⁠) are numerically determined (same data for bullets of Fig. 14); red line is $n = 1$ [i.e., $β_{F} (E_{m}) = β_{T}$ ]. The field regime is blue region; the thermal regime is orange region; the thermal-field regime is where $| n - 1 | \leq 0.05$ ⁠.

FIG. 16.

View large Download slide

Integrand of $N (n, s)$ of Eq. (96) normalized to unity for each value of $n$ ⁠. Yellow dots show the location of the maximum. White lines at $z = 0$ correspond to $E = μ$ ⁠, and at $z = s$ to $E = μ + ϕ$ ⁠. The thermal-field regime is mostly between the white lines; to the left is the field regime, and to the right is the thermal regime.

The behavior of the (normalized) integrands for the exact $d J (E)$ and its approximation by the FN and quadratic approximations to $θ (E)$ in Fig. 17, as well as the comparison to the integrand when the $L θ$ approximation of Eq. (73) is used as shown in Fig. 18, gives the mistaken impression that the approximations are equivalent in terms of their estimates of the integrated current density $J (F, T)$ ⁠: the curves in those figures are normalized because the maximum of the current integrands varied widely for the FN and quadratic approximations (the maximum of the $L θ$ approximation is by definition the same as the exact approach), and so viewing normalized integrands aids in understanding their evolution with both $F$ and $T$ but obscures stark differences in their overall magnitudes. Tables IV and V show the maximum and the ratio of the integrated functions, respectively: the latter ratios give a measure of the error associated with the $L θ$ approximation in calculating $J (F, T)$ ⁠. When the differences in integrand maximum height are understood to be generally large, then it is clear that the $L θ$ substantially outperforms either FN or Quad in the TF region, and still performs better than either in field or thermal regimes, respectively.

FIG. 17.

View large Download slide

Current density integrand $d J (F, T)$ normalized to its maximum for the SN barrier as a function of $E$ for $F = 2$ eV/nm, $μ = 7$ eV, and $Φ = 4.5$ eV. Thermal regime is shaded orange, field regime is shaded blue, with thermal-field regime the white space between. Integrand for the exact method is filled gray. FN is Eq. (87), and Quad is Eq. (83).

FIG. 18.

View large Download slide

Solid line is the exact integrand $d J (E)$ of (87). Dashed line is the $L θ$ approximation of Eq. (73). Parameters are $μ = 7$ eV, $Φ = 4.5$ eV, and $F = 2$ eV/nm. Legend is the value of $T$ in (K).

TABLE IV.

Ratio of the maximum of the integrand dJ_x(E) to dJ(E_m), where x ∈ (Lθ, FN, Quad), for the T of Figs. 17 and 18 for μ = 7 eV, Φ = 4.5 eV, and F = 2 eV/nm.

T [K]	Lθ	FN	Quadratic
673	1.0000	0.9008	338.1524
1173	1.0000	7.1085	93.6029
1673	1.0001	14.7393	1.0521
2173	1.0000	9.3149	0.9747
2673	1.0000	6.6280	0.9652

T [K]	Lθ	FN	Quadratic
673	1.0000	0.9008	338.1524
1173	1.0000	7.1085	93.6029
1673	1.0001	14.7393	1.0521
2173	1.0000	9.3149	0.9747
2673	1.0000	6.6280	0.9652

TABLE V.

Ratio of $J_{x} = \int d J_{x} (E)$ to $J = \int d J (E)$ ⁠, where x ∈ (Lθ, FN, Quad), for the temperatures of Figs. 17 and 18 for μ = 7 eV, Φ = 4.5 eV, and F = 2 eV/nm.

T (K)	Lθ	FN	Quadratic
673	1.0305	0.9277	375.6567
1173	1.3969	3.9578	41.7142
1673	1.3257	8.2602	2.0364
2173	1.1009	6.4468	1.1889
2673	1.0574	4.9732	1.0733

T (K)	Lθ	FN	Quadratic
673	1.0305	0.9277	375.6567
1173	1.3969	3.9578	41.7142
1673	1.3257	8.2602	2.0364
2173	1.1009	6.4468	1.1889
2673	1.0574	4.9732	1.0733

2. Numerical evaluation

In the original GTF (⁠ $o$ GTF) approach,^5,8,103 the non-linearity of $θ (E)$ was approximated by considering a cubic polynomial with coefficients determined by analytically well-understood values of $θ (μ), θ (μ + ϕ) = 0, β_{F} (μ)$ ⁠, and $β_{F} (μ + ϕ)$ using conventional FN and quadratic barrier methods. Its numerical implementation^5,8 proceeded rapidly as follows for variations in temperature:

For a given $F$ ⁠, define the temperatures $T_{min}$ and $T_{max}$ as the smallest and largest temperatures for which $n (F, T) = 1$ ⁠.
For $T < T_{min}$ ⁠, $E_{m} = μ$ ⁠; for $T > T_{max}$ ⁠, $E_{m} = μ + ϕ$ ⁠; for $T_{min} < T < T_{max}$ ⁠, solve $β_{F} (E_{m}) = β_{T}$ ⁠.
Evaluate $n (F, T) = β_{T} / β_{F} (E_{m})$ and $s (F, T) = θ (E_{m}) + β_{F} (E_{m})$ $[E_{m} - μ]$ ⁠.
Approximate $N (n, s)$ using Eq. (97).
Evaluate $J (F, T)$ using Eq. (96).

The same algorithm is used for variations in the field, although the determination of

F_{max}

and

F_{min}

from Eq. requires numerical rather than analytic methods. The

o

GTF approach is computationally rapid with low numerical overhead and reproduces the canonical FN and RLD equations in the appropriate limits, as it was by construction designed to do with its choices of

E_{m}

⁠. However, as shown by Kyritsakis and Xanthakis,¹¹⁹ the approximation suffers because the approximations to

E_{m}

depart from those numerically found particularly in or near the TF regime, causing small cusps apparent in the transitions between thermal to thermal-field and thermal-field to field regimes (see Fig. 5 of Ref. 119). Even though

θ (E)

is only weakly non-linear, current density depends on exponentials of it, and so not finding

E_{m}

exactly coupled with the weak non-linearity of

θ

and the small variations in

β_{F} (E)

has visible consequences. Correcting those cusps becomes desirable, even though the

o

GTF method already outperforms the canonical equations in terms of accuracy in the TF regime but recovers them in either the limiting cases of the thermal or field regimes.

A reformulated (⁠ $r$ GTF) approach^92,120 was subsequently developed that provided a method for finding $E_{m}$ directly, from which $θ (E_{m})$ and $β_{F} (E_{m})$ were then calculated. Doing so results in superior estimates of $n (F, T)$ and $s (F, T)$ ⁠. The present work is based $r$ GTF but with further modifications to improve computational expediency, specifically by insisting that $E_{m}$ ⁠, $θ (E_{m})$ ⁠, and $β_{F} (E_{m})$ all be found in a computationally rapid manner. An impediment to rapidity is the complexity of Eqs. (68) and (84) and a requirement for their evaluation for $y > 1$ (e.g., $E_{m} > μ + ϕ$ ⁠) in the pure thermal limit. The improvement makes use of both $θ$ and $β_{F}$ being relatively smooth functions in both $E (y)$ and $y$ ⁠. A computationally expedient method is drawn from an analysis of $θ (E)$ ⁸⁶ for various barriers. The upgrades improve on the accuracy of $r$ GTF, with an ancillary benefit that straightforward extensions to barriers beyond the SN barrier, e.g., MIM, depletion, prolate spheroidal, and space–charge modified barriers, become possible. The $r$ GTF enhanced method is computationally implemented as follows:

Determine $θ (y)$ from Eq. (68), and $E (y) = μ + Φ - \sqrt{4 Q F} / y$ for four equispaced values in $E (y)$ between $y_{m i n} = \sqrt{4 Q F} / (μ$ $+ ϕ)$ and $y_{m a x} = 1$ ⁠, where $ϕ = Φ - \sqrt{4 Q F}$ ⁠.
From $θ (y)$ and $E (y)$ ⁠, fit $θ [y (E)]$ by a cubic equation in $(μ + ϕ - E)$ ⁠, where $V_{o} \equiv μ + ϕ$ ⁠, which will generate three $C_{j}$ (⁠ $C_{0}$ is by definition zero).
Find $θ (E_{m})$ and $β_{F} (E_{m})$ for Steps 4–6 below by

$\begin{matrix} θ (E_{m}) \Leftarrow \sum_{j = 1}^{3} C_{j} | μ + ϕ - E_{m} |^{j}, \\ β_{F} (E_{m}) \Leftarrow \sum_{j = 1}^{3} j C_{j} | μ + ϕ - E_{m} |^{j - 1} . \end{matrix}$
(100)
Take $θ (μ + ϕ - E) = - θ (| μ + ϕ - E |)$ (antisymmetric) when $E > μ + ϕ$ ⁠.
Approximate $N (n, s)$ using Eq. (97).
Evaluate $J (F, T)$ using Eq. (96).

There is no inconsistency in the

L θ

reliance on

θ (E)

being approximated locally by a linear function in

E

and the

r

GTF approach approximating

θ

as a cubic in

(μ + ϕ - E)

⁠: both

n (F, T)

and

s (F, T)

depend only on the value and first derivative of

θ (E)

at

E = E_{m}

in order to approximate

N (n, s)

⁠. Nevertheless, outside the thermal-field region, the

r

GTF method replicates the FN (⁠

n \to \infty

⁠) and RLD (⁠

n \to 0

⁠) behavior in keeping with Eq. and as directly indicated in Fig. 19. In the absence of the

Σ (n)

and

Σ (1 / n)

factors in Eq. ,

J_{M G} (F, 0)

and

J_{R T D} (F, T)

are reproduced, but it is seen that the neglect of the

Σ

-factors is consequential by the separation between the

°

-data points and the (⁠

⋄, □

⁠)-data points.

FIG. 19.

View large Download slide

Comparison of $N (n, s)$ (green $°$ ⁠) to $e^{- n s}$ (red $□$ ⁠) and $n^{2} e^{- s}$ (blue $⋄$ ⁠) to show consequences of neglecting $Σ (n)$ and $Σ (1 / n)$ in Eq. (98). To show convergence of the latter two to the canonical equations, $N_{f n} = J_{M G} / A_{R L D} T^{2}$ and $N_{r l d} = J_{R L D} / A_{R L D} T^{2}$ are shown as dashed lines (cyan and orange, respectively).

Although the

r

GTF method substantially improves agreement with numerical evaluations of

J (F, T)

in the thermal-field regime (⁠

n \approx 1

⁠) compared to

o

GTF, nevertheless, as suggested by Figs. 20 and 21, the

r

GTF method overestimates the exact

J (F, T)

⁠. Such accuracy may be sufficient: recall that the FN equation suffers comparable or worse departures from the numerical evaluation of Eq. but is nevertheless widely used because it is computationally rugged. If greater accuracy is required but numerical integration is too computationally costly, then a correction developed for

r

GTF can be utilized: a Lorentz factor¹²⁰ reduces the discrepancy at the peaks of Fig. 21. The method involves using a small number of evaluations to approximate the ratio of the

L θ

evaluation to the exact approach (that is,

J_{n s} / J

⁠) as a function of, say, temperature, by a sum of Lorentian terms, as shown in Fig. 21, where

\begin{matrix} \frac{d J_{n s}}{d J} \approx 1 + h_{1} L_{1} (T / T_{o}) + h_{2} L_{2} (T / T_{o}), \\ L_{j} (x) \equiv \frac{1}{1 + p_{j} {(x - ϵ_{j})}^{2}}, \end{matrix}

(101)

as shown in Fig. 21, where

J = \int d J \approx \sum_{j} d J_{j} Δ E

⁠. The parameters are conveniently extracted by finding the location of the maxima of the two peaks and using ad hoc values of

p_{j}

to find, and then invert the equation

[\begin{matrix} 1 & L_{2} (ϵ_{1}) \\ L_{1} (ϵ_{2}) & 1 \end{matrix}] [\begin{matrix} h_{1} \\ h_{2} \end{matrix}] = [\begin{matrix} H (ϵ_{1}) \\ H (ϵ_{2}) \end{matrix}],

(102)

where

L_{j} (ϵ_{j}) \equiv 1

and

H (ϵ_{j}) = [J_{n s} (F_{j}, T_{j}) / J (F_{j}, T_{j})] - 1

⁠. As seen in Fig. 21, even crude and unoptimized Lorentzians provide reasonably accurate corrections. Although agreement can be improved by an optimization algorithm with the addition of a window function,¹²⁰ such an approach shall not be further considered below because

total current is often represented on log-plots for which the changes are only weakly discernible,
the computational overhead to find the best Lorentzians for each $F$ may be undesirable and, therefore, ill-suited for simulations demanding millions of evaluations of $J (F, T)$ in a beam simulation code, and
the $r$ GTF method using $E_{m}$ but not the Lorentzian correction is already a better approximation than either the FN or RLD canonical equations, both of which are employed without hesitation in the literature already, but which depart significantly from $J (F, T)$ [Eq. (72)] by comparison.

FIG. 20.

View large Download slide

A contour plot representation of Fig. 18 showing the actual width of the $d J (E)$ integrand (filled contour lines) compared to its $L θ$ approximation (dashed lines of the same color). The width of the integrand in the thermal-field region between the vertical white lines at $E = μ$ and $E = μ + ϕ$ shows why the $L θ$ approximation results in a larger estimate of $J = \int d J$ ⁠.

FIG. 21.

View large Download slide

The ratio $J_{s n} / J$ as a function of $T / T_{o}$ (blue $∙$ ⁠) and $n (F, T)$ (red $□$ ⁠) for the $x$ axis. “L-fit” refers to Eq. (101) with $(h_{1}, p_{1}, a_{1})$ $= (0.3808, 1000, 0.8422)$ ⁠, $(h_{2}, p_{2}, a_{2}) = (0.6634, 125, 1.125 87)$ ⁠, and $T_{o}$ $= 1250$ K. Symbols (⁠ $∙, □$ ⁠) occur at the same $J_{s n} / J$ in each $x$ -coordinate.

Even in the absence of a Lorentzian correction, the $r$ GTF method provides a good account of experimental data. A quintessential example is afforded by considering the measurements of Dyke et al.,³⁴ who measured TF emission from a tungsten needle. The emitter was subjected to six different elevated temperatures and the field emission current measured. Comparing those values to the current drawn from a room temperature tungsten needle resulted in the data points shown in Fig. 22. Using the fields provided by Dyke et al.³⁴ as given (unmodified) gives an equivalent behavior, but the inference of the fields at the apex of the tungsten emitters was based on theoretical and geometrical models used by those authors. It is, therefore, to be expected that the actual fields may depart in unknown ways from the inferred fields due to local field enhancement, work function variations over the apex of the emitter, departure of the idealized emitter from the actual emitter, the integrated current density over the emitter and its relation to current density at the apex,¹²⁰ and other probable complications. On the assumption that the apex field on the emitter bears a reasonably linear relationship to the anode voltage in an experimental configuration, assume that the apex field departs slightly from the values reported by Dyke et al., and that the actual fields are the form $F = a F_{d b} + b$ ⁠, where $F_{d b}$ are the values obtained by Dyke et al.³⁴ The slope/intercept values $(a, b)$ are found by adjusting $F$ in $J (F, T)$ for each data set until the theoretical curve matches the experimental data: when the fitted $F$ are plotted as a function of $F_{d b}$ ⁠, the relationship is to a good approximation linear, with $F$ $= 1.2559 F_{d b} - 0.5027$ in units of [eV/nm], as used in Fig. 22, for which the highest temperature cases for each line are approaching $n \approx 1$ from above. Evaluation of $n (F, T)$ directly shows that at the highest currents for each data set, the transition to the TF regime has occurred. As a result, it is demonstrated that usage of $J_{M G} (F, 0)$ alone can be in error by an order of magnitude when temperature effects are neglected, apart from the neglect of thermal contributions from other parts of the needle away from the apex.

FIG. 22.

View large Download slide

Comparison of experimental thermally assisted field emission from a tungsten needle (symbols, digitally extracted from Fig. 22 of Ref. 34) where the fields $F$ used in the $r$ GTF lines are related to the fields $F_{b d}$ of Barbour et al. by $F = 1.2559 F_{d b} - 0.5027$ ⁠. The data are labeled by $F_{b d}$ ⁠.

D. Photoemission

1. Extended Fowler–DuBridge model

A straightforward model of photoemission was developed by Fowler¹²¹ and modified by DuBridge¹²² by making use of the framework familiar from the field emission derivation. In it, a photoexcited electron was assumed to absorb a photon energy of

ℏ ω

such that the entire energy of the photon contributed to the forward energy into the surface barrier. Electrons for which

E

+ ℏ ω > μ + Φ

were then emitted. The approximation is reasonable for electrons emitted at the threshold. In examining the derivation of the FN equation, it is seen that the approximation

E \to E + ℏ ω

is tantamount to

s (F, T)

changing from a positive quantity to a negative one, and so for photoemission,

N (n, s) \to N (n, - | s |)

⁠. It is, therefore, more convenient for exposition to redefine

s

in the case of photoemission to make it a positive quantity and thereby speak of

N (n, - s)

as applying to photoemission by assuming

the initial state of the photoexcited electron is at the Fermi level $μ$ ⁠,
the emission distribution is peaked near $E_{m} \approx μ + ϕ$ ⁠, and so
take $- s = β_{F} [μ + ϕ - (μ + ℏ ω)] = β_{F} (ℏ ω - ϕ)$ ⁠.

An evaluation of the sum of

N (n, - s)

with

N (n, s)

by grouping terms in their series expansions¹⁰³ shows that to a good approximation,

N (n, s) + N (n, - s) = \frac{1}{2} {(n s)}^{2} + ζ (2) (n^{2} + 1) .

(103)

Because

N (n, s)

is exponentially small,

N (n, - s)

is mostly given by the right hand side. As quantum efficiency is proportional to emitted current, therefore

Q E \propto T^{2} N (n, - s)

⁠, and therefore, because

n = β_{T} / β_{F}

⁠,

k_{B} T = 1 / β_{T}

⁠, and

n s = β_{T} (ℏ ω - ϕ)

it follows

Q E \propto \frac{1}{2} {(ℏ ω - ϕ)}^{2} + ζ (2) (β_{F}^{- 2} + β_{T}^{- 2}),

(104)

where

ζ (2) = π^{2} / 6

⁠. Measurements of photoemission near the threshold give evidence of the behavior of Eq. , although in the laboratory, fields on the photoemitter tend to be so small (10–100 kV/m) that

1 / β_{F}

is often negligible (even in accelerators where fields are on the order of 1–100 MV/m corresponding to

F = 0.001 - 0.1

eV/nm,

β_{F}

may be large) and allow for the neglect of

1 / β_{F}^{2}

in Eq. when examining measured data: as a comparison, for copper-like conditions and using the parabolic barrier approximation to the SN barrier of Eq. gives

1 / β_{F}^{2} = 6.433 94 \times 10^{- 6}

eV

^{2}

for

F = 0.01

eV/nm, compared to

1 / β_{T}^{2} = 6.683 27 \times 10^{- 4}

eV

^{2}

for

T = 300

K, or about two orders of magnitude larger. Copper is a common and well-characterized photocathode material^14,123 for which QE data are available,¹²⁴ as shown in Fig. 23. To make a comparison to theory, the theory curves are generated by

Q E_{t h} (ω) = Q E_{\exp} (ω_{o}) [\frac{{(ℏ ω - ϕ)}^{2} + 2 ζ (2) k_{B} T^{2}}{{(ℏ ω_{o} - ϕ)}^{2} + 2 ζ (2) k_{B} T^{2}}],

(105)

for

T = 300

K and

ℏ ω = 5

eV, for which

2 ζ (2) {(k_{B} T)}^{2} = 2.198 71

\times 10^{- 3}

eV

^{2}

⁠. It is evident that the model performs well near the threshold for a particular value of

ϕ = 4.36

eV, but it is also evident that departures occur at higher energy, the causes of which shall be considered below.

FIG. 23.

View large Download slide

The quantum efficiency of a cleaned copper surface. Circles (⁠ $°$ ⁠) correspond to experimental measurements (data courtesy of D. Dowell, SLAC). Lines correspond to Eq. (104) neglecting $1 / β_{F}^{2}$ and assuming $T = 300$ K and $ℏ ω_{o} \approx 5$ eV for different values of $ϕ$ shown in the legend (in eV). The termination of the theory lines occurs where $ℏ ω \approx ϕ$ ⁠.

The same analysis applied to high quantum efficiency semiconductors like Cs $_{3}$ Sb show deviations from the metal model, as in Fig. 24 for the data of Spicer.¹¹ That it does so reflects that for electron energies near the barrier maximum, Eq. (72) continues to describe emission, but subject to modifications that must be accounted for, of which there are several:

semiconductors are characterized by an electron affinity $E_{a}$ (height of the barrier above the conduction band minimum) and a band gap $E_{g}$ ⁠;
photoexcited electrons are not in thermal equilibrium (Fermi–Dirac distribution) so that scattering events between them and other electrons or the lattice cause energy relaxation that renders some (in the case of semiconductors) or most (in the case of metals) electrons no longer eligible for emission;¹²⁵ and
the three-dimensional nature of the Fermi distribution that was buried in $f (k_{x})$ of Eq. (42) must be taken into account.

A primary difference between semiconductors and metals is the presence of a band gap for the former, which complicates the scattering and changes the basic three-dimensional model.¹²⁶ A computationally convenient account, though, profits from considering the limit of

β_{T} \to \infty

(low temperature).

FIG. 24.

The quantum efficiency of Cs 3Sb ( °) digitally extracted from Fig. 8 of Spicer.11. Lines correspond to Eq. (104) neglecting 1 / β F 2 and assuming T = 300 K and ℏ ω o ≈ 2.5 eV for different values of ϕ shown in the legend (in eV). The termination of the theory lines occurs where ℏ ω ≈ ϕ.

View large Download slide

The quantum efficiency of Cs $_{3}$ Sb (⁠ $°$ ⁠) digitally extracted from Fig. 8 of Spicer.¹¹. Lines correspond to Eq. (104) neglecting $1 / β_{F}^{2}$ and assuming $T = 300$ K and $ℏ ω_{o} \approx 2.5$ eV for different values of $ϕ$ shown in the legend (in eV). The termination of the theory lines occurs where $ℏ ω \approx ϕ$ ⁠.

With elegant simplicity, the FD model of Eq. captures the behavior of quantum efficiency near the threshold, but the approach additionally benefits from being able to describe modifications brought on by changes in the photoexcited electron population. As an example, Rashba spin–orbit coupling causes photoexcited electrons to separate into different spin chirality populations and provides a path to spin-polarized currents different from the usage of GaAs (currently the only technology practical for producing high intensity spin-polarized electron beams for accelerators capable of meeting the needs of nuclear and particle physics experiments¹²⁷). An analysis replacing the standard parabolic electron dispersion with one characterizing a Rashba spintronic material (RSM) results in a modified FD model (using methods recognizable from those herein) to find¹²⁸

Q E \propto \frac{1}{2} {(ℏ ω - ϕ)}^{2} + E_{R} (ℏ ω - ϕ) - \frac{E_{R}^{3}}{6},

(106)

where

E_{R}

is a factor affecting the apparent work function through the Rashba spin–orbit coupling (RSOC) effect and generally a fraction of an eV.

2. Moments model

A three-step model of emission, in which the processes of absorption, transport, and emission are modeled by various components, is a useful method to describe both photoemission and secondary emission processes.^{8,12,13,16,123} For photoemission, the steps are

Absorption (⁠ $A$ ⁠), in which incident light (typically a laser in photoinjectors) illuminates the surface of a photocathode, where some of the light is reflected (governed by $R (ω)$ ⁠) and some absorbed but attenuated as it passes into the photocathode material [governed by $δ (ω)$ ];
Transport (⁠ $T$ ⁠), in which a photo-excited electron absorbs the photon energy $ℏ ω$ ⁠, and the excited electron moves in a random direction (specified by a polar angle $θ$ ⁠) which may take it back to the surface, during which time it may suffer scattering events (governed by a scattering rate $τ$ ⁠); and
Emission (⁠ $E$ ⁠), in which the photo-excited electron escapes (which is governed by a transmission probability $D (E)$ ⁠).

The simple moments model¹²⁶ adopts that sequence of events but uses a momentum-based rather than energy-based framework to evaluate an integral for the current density of the form

J (ℏ ω) = \frac{2}{{(2 π)}^{3}} ∭ (\frac{ℏ k_{z}}{m}) A \times T \times E d k_{z} d {\vec{k}}_{⊥} .

(107)

The quantum efficiency is then the ratio of

J (ℏ ω)

when all the impediments to emission are in play to a current density for which all photoexcited electrons suffer no impediment to emission (infinite mean free path and unitary transmission probability). The coefficient of 2 that for thermal and field emission accounted for electron spin is unnecessary here, but included because taking ratios eliminates it anyway, and to retain a connection to the other emission equation formulations. The mathematical framework can be reconfigured as moments of a distribution function

f (\vec{k})

where the various processes depend on position, for which the moments are

M_{n} (k_{j}) \equiv \frac{2}{{(2 π)}^{3}} \int k_{j}^{n} f (\vec{k}) d^{3} k,

(108)

where

d^{3} k = 2 π k^{2} \sin θ d θ

⁠,

k_{j}

is a Cartesian component of

\vec{k}

⁠, e.g.,

k_{z} = k \cos θ

⁠, and a factor of 2 in the numerator of the coefficient accounts for electron spin. The distribution implicitly accounts for the electrons that successfully escape the surface. Quantum efficiency then becomes

Q E \equiv [1 - R (ω)] \frac{M_{1} (k_{z})}{2 {M_{1} (k_{z}) |}_{D = 1, f_{λ} = 1}},

(109)

where

R (ω)

accounts for losses due to reflected photons,

D (k_{z})

is the transmission probability as a function of forward momentum

k_{z}

⁠,

f_{λ}

is a factor that accounts for losses due to scattering events, and the factor of 2 in the denominator accounts for the fact that half of the photoexcited electrons on average are traveling away from the surface rather than toward it. In the moments model, the number density is the 0th moment, and the current density is the product of the 1st moment of

k_{z}

with

(ℏ / m)

⁠. The same equation is seen to be related to Eq. by taking

ℏ ω \to 0

⁠, but a substantial modification arises for photoemission: the probability that an electron has a total energy

E

before the absorption of a photon is still governed by the Fermi–Dirac distribution

f_{F D} (E)

⁠, but now the probability that its final state is free must be taken into account by a factor of

[1 - f_{F D} (E + ℏ ω)]

⁠, for which the low temperature limit results in the introduction of Heaviside step functions

Θ (x)

(0 for negative argument, 1 for positive argument),

f_{F D} (E) [1 - f_{F D} (E + ℏ ω)] = Θ (μ - E) Θ (E - μ + ℏ ω),

a restriction that places limits on the momentum integrations. For low fields, the transmission probability

D (k_{z})

prevents electrons escaping for which

ℏ^{2} k_{z}^{2} / 2 m < μ + ϕ

for metals with a Fermi level

μ

and an effective work function

ϕ = Φ - \sqrt{4 Q F}

⁠, with a similar equation for semiconductors with a positive electron affinity

E_{a}

and band gap

E_{g}

⁠. Letting

ℏ k_{z} / m = (ℏ k / m) \cos θ

⁠, then a maximum angle

θ_{m}

exists for which emission is allowed and is given by

(E

+ ℏ ω) \cos^{2} θ_{m} = μ + ϕ

for metals, or

E \cos^{2} θ_{m} = E_{a}

for semiconductors,

{\cos θ_{m} |}_{m e t a l} = {[\frac{μ + ϕ}{E + ℏ ω}]}^{1 / 2},

(110)

{\cos θ_{m} |}_{s e m i} = {[\frac{E_{a}}{E}]}^{1 / 2},

(111)

but with the caveat that energy is being treated differently for metals and semiconductors: in Eq. ,

E

is the energy of the electron before excitation, whereas in Eq. , it is more convenient to let

E

be the energy after excitation, a distinction required by how the integrals behind

Q E

are set up in SSM. Consequently, although the metal and semiconductor versions of the various equations appear to use the same

E

⁠, they do not, and how the calculation proceeds must take that into account.¹²⁹ To reinforce that distinction,

\tilde{E}

will designate the energy of the photoexcited electron, and so

\tilde{E} = E + ℏ ω

for metals, or

\tilde{E} = E

for semiconductors in what follows.

Working in polar coordinates then, with total energy

E

= ℏ^{2} k^{2} / 2 m

and

d E = (ℏ^{2} k / m) d k

allows the SMM to be recast in terms of

E

and

x_{m} = \cos θ_{m}

as

Q E_{m e t a l} = \frac{(1 - R (ω))}{ℏ ω (2 μ - ℏ ω)} \int_{Δ}^{μ} E d E F_{λ} (x_{m}, E),

(112)

Q E_{s e m i} = \frac{(1 - R (ω))}{{(ℏ ω - E_{g})}^{2}} \int_{E_{a}}^{Δ_{s}} E d E F_{λ} (x_{m}, E),

(113)

F_{λ} (x_{m}, E) \equiv \int_{x_{m}}^{1} D (E s^{2}) f_{λ} (s, E) s d s,

(114)

where

Δ_{m} = μ + ϕ - ℏ ω

for metals and

Δ_{s} = ℏ ω - E_{g}

for semiconductors,

x_{m} = \cos θ_{m}

⁠, and the denominators of the coefficients are the numerator integrals configured such that all photoexcited electrons escape. For metals and semiconductors, respectively, the denominators go as

\begin{aligned} 2 \int_{μ - ℏ ω}^{μ} E d E \int_{0}^{π / 2} \sin θ d θ = ℏ ω (2 μ - ℏ ω), \\ 2 \int_{0}^{ℏ ω - E_{g}} E d E \int_{0}^{π / 2} \sin θ d θ = {(ℏ ω - E_{g})}^{2}, \end{aligned}

with

ϕ \to 0

and

F_{λ} \to 1

(corresponding to all photoexcited electrons escaping). The scattering factor

f_{λ}

is such that an electron photoexcited a distance

z

from the surface will with equal probability be traveling in any direction from its excitation point, and so the distance it must travel to get back to the surface in the direction of

\vec{k}

is

z / \cos θ

for

θ < π / 2

(larger

θ

are traveling away from the surface and so do not contribute to QE). Let the mean free path between scattering events be

ℓ (\tilde{E})

⁠, where

\tilde{E}

is the energy of the photoexcited electron,

ℓ = v (\tilde{E}) τ = (τ / m) \sqrt{2 m \tilde{E}}

and

τ (\tilde{E})

is the relaxation time describing scattering and is energy-dependent.

First, consider metal photoemission conditions. Let

Δ_{m} \to Δ = μ + ϕ - ℏ ω

⁠, and replace

\tilde{E} = E + ℏ ω

⁠, where

E

is the energy of the electron before photoexcitation. Further, approximate

D (E s^{2})

as a step function, so that it is simply unity over the range of

s

-integration. Electron–electron scattering in metals reapportions the initial energy of photoexcited electron among the daughter electrons such that to a good approximation, a single scattering event is enough to render both daughter electrons below the surface barrier: the fatal approximation, therefore, assumes that any scattering event is “fatal” to emission, although photons of sufficiently high energy may allow a single scattering event to be survivable for metals.^125,130 If the absorption of photons in a material is governed by a laser penetration depth

δ (ω)

such that the probability of absorption at a distance

z

from the surface is proportional to

\exp (- z / δ)

⁠, then the fraction that makes the journey to the surface is given by

\begin{aligned} f_{λ} (\cos θ, E) & \equiv \frac{1}{δ} \int_{0}^{\infty} \exp [- \frac{z}{δ} (1 + \frac{p}{\cos θ})] d z \\ = \frac{\cos θ}{\cos θ + p (E)}, \end{aligned}

(115)

where

p (E) = δ (ω) / ℓ (\tilde{E})

⁠. The integral in Eq. is then analytic and given by

F_{λ} (x, p) = \frac{1}{2} (1 - x) (1 + x - 2 p) - p^{2} \ln (\frac{x + p}{1 + p}),

(116)

where the subscript on

x_{m}

and the energy arguments of

x_{m}

and

p

are momentarily suppressed. Incorporation of

F_{λ} (x, p)

into the integrals of Eqs. and is non-trivial by virtue of the energy dependence of both

x_{m} (E)

and

p (E)

⁠, and, therefore, numerical means are required if Eq. is used, but a leading order approximation with minor computational overhead is possible. First, approximate

p (E) \approx p (μ)

throughout (electrons photoexcited from closer to the Fermi level will tend to dominate). Second, observe that

F_{λ} [x_{m} (Δ), p] = 0

exactly because

x_{m} (Δ) = 1

(the lowest energy photoexcited electron to escape is normal to the surface and has no transverse energy components), and so the roughly linear behavior of

F_{λ}

may be approximated by

F_{λ} [x_{m} (E), p (E)] \approx F_{λ} [x_{m} (μ), p (μ)] [\frac{E - Δ}{μ - Δ}],

(117)

in Eq. , resulting in the computationally expedient replacement,

Q E_{m e t a l} \approx \frac{N_{ω}}{6} (ℏ ω - ϕ) (3 μ + ϕ - ℏ ω) F_{λ} [x_{m} (μ), p (μ)],

(118)

where

N_{ω} = [1 - R (ω)] / [ℏ ω (2 μ - ℏ ω)]

⁠. The connection to the Fowler–DuBridge equation is instructive. It is obtained by first letting

f_{λ} (\cos θ, E) \approx \cos θ / [1 + p (μ)]

⁠, for which

F_{λ} [x_{m} (E), p (μ)

\to [1 - x_{m} {(E)}^{3}] / [3 (1 + p (μ))]

⁠. Using that approximation in the above integrals and expanding

x_{m} (μ)

to first order in

(ℏ ω - ϕ)

results in

F_{λ} [x_{m} (μ), p] \int_{Δ}^{μ} s d s \approx μ {(ℏ ω - ϕ)}^{2} / [2 (μ + ϕ) (1 + p)]

⁠, which, when inserted into Eq. , results in [compare Eq. of Ref. 126]

Q E_{m e t a l} \sim \frac{1 - R (ω)}{4 [1 + p (μ)]} \frac{μ {(ℏ ω - ϕ)}^{2}}{ℏ ω (2 μ - ℏ ω) (μ + ϕ)},

(119)

and is recognized as reproducing the dependence of the Fowler–DuBridge equation on

{(ℏ ω - ϕ)}^{2}

to leading order. A comparison of the exact evaluation of Eq. with the computationally expedient approximation of Eq. and the near threshold approximation of Eq. is shown in Fig. 25 compared to the data of Dowell et al. from Fig. 23. The relaxation time

τ (E)

and laser penetration depth

δ (ω)

were implicitly held fixed, but both change as a consequence of photon energy

ℏ ω

in a manner that may allow better correspondence away from threshold because

p (μ) = δ (ω) / l (\tilde{μ}) = δ (ω) / l (μ + ℏ ω)

as per Eq. .

FIG. 25.

View large Download slide

The measured photoemission data of Fig. 23 for copper (data courtesy of D. Dowell, SLAC) compared to Eq. (118) (red, “theory”) and Eq. (119) (blue, “FD limit”). The normalization factor for the theory lines is found by demanding that $Q E (ℏ ω_{o}) |_{t h e o r y} = Q E (ℏ ω_{o}) |_{e x p} = 0.000 135$ for $ℏ ω_{o} = 4.88$ eV and $ϕ = 4.34$ eV.

Second, consider semiconductor photoemission conditions. Again, let

Δ_{s} \to Δ = ℏ ω - E_{g}

⁠,

x_{m} = \sqrt{E_{a} / E}

⁠, and use the trapezoidal approximation to approximate the integral behind

F_{λ} (x_{m}, E)

to give

\begin{aligned} F_{λ} (x, E) & = \int_{x_{m}}^{1} \frac{x^{2}}{x + p} D (E x^{2}) d x \\ \approx \frac{2}{1 + p} \frac{{(x - E_{a})}^{5 / 2}}{(\sqrt{x} + \sqrt{E_{a}}) {(\sqrt{x} + \sqrt{x - E_{a}})}^{2}} . \end{aligned}

(120)

Likewise evaluating Eq. for

Q E_{s e m i}

with Eq. and using a trapezoidal approximation again gives

\begin{matrix} Q E_{s e m i} \approx [\frac{1 - R (ω)}{1 + p}] h (s), \\ h (s) \equiv \frac{s^{5}}{u^{2} (1 + u) {(s + u)}^{2}}, \end{matrix}

(121)

where

s \equiv \sqrt{(ℏ ω - E_{a} - E_{g}) / E_{a}}

⁠,

u \equiv \sqrt{1 + s^{2}}

⁠,

p (E)

is evaluated at

E = E_{a}

⁠, and the form here corrects a typographical error in prior accounts^8,32,126 that inadvertently neglected the power of

2

on

{(s + u)}^{2}

in the denominator of

Q E_{s e m i}

(the numerical calculations of

Q E_{s e m i}

in the figures and code of those references used the correct formula). Observe that asymptotically as

ℏ ω \to \infty

⁠, then the

s

and

u

dependent fraction in Eq. approaches unity. Equation is to be contrasted with the parametric form suggested by Spicer¹¹ that can be recast using the present notation as¹³¹

Q E_{s p i c e r} = G (ℏ ω) (\frac{s^{3}}{s^{3} + γ}),

(122)

where

G

corresponds to Spicer’s coefficient

G (h ν)

⁠: because

h (s) \to 1 / 4

⁠, as

s \to \infty

⁠,

G

corresponds to a quarter of the coefficient in square brackets of Eq. . For a more modern (and complex) version, see Ref. 15 and references therein. Spicer’s parametric representation is designed for threshold emission conditions, as seen in the comparison of

h (s)

with Spicer’s

h_{s p} (s) = G_{o} s^{3} / [4 (s^{3} + γ)]

form in Fig. 26.

FIG. 26.

View large Download slide

Comparison of $h (s)$ [Eq. (121)] vs the parametric approximation due to Spicer [Eq. (122)] for $G_{o} / 4 = h (4)$ and $γ$ as shown in the legend.

A comparison of Eq. (121) to the photoemission data of Spicer,¹¹ as well as Taft and Philipp¹³² is shown for the presumed parameters of $E_{g} = 1.6$ and $E_{a} = 0.4$ eV: a single set of parameters describes both experimental data sets well, even though the cited sources reported differing values for those parameters compared to each other. As with the analysis for copper data, the relaxation time $τ$ and reflectivity $R (ω)$ are held fixed: models exist for their variation¹²⁶ that are expected to improve the correspondence, but as a simple calculational tool, the SSM model performs quite well and has the advantage that once $R (ω)$ and $p (E_{a})$ are fixed by other means, only one photoemission data point is then required to obtain $Q E_{s e m i}$ over a range of frequencies for a given $(E_{g}, E_{a})$ set of values, in contrast to $Q E_{s p i c e r}$ [Eq. (122)], which has both $G (ℏ ω)$ and $γ$ as connected fitting factors that differ for each experimental data set and suggest differing values of $E_{a}$ and $E_{g}$ as well as $γ$ and $G_{o}$ ⁠.³² In contrast, the SMM approach separates the determination of $E_{g}$ and $E_{a}$ [which, if not input parameters, are determined by the shape of $h (s)$ ] from the separate determination of $R (ω)$ and $p (E)$ for which the ratio $(1 - R) / (1 + p)$ can be set, separately from $h (s)$ ⁠, using a single experimental value of $Q E$ ⁠. For example, the red line in Fig. 27 is $h (s)$ ⁠, and the green dot corresponds to $Q E_{e x p} (ω_{o}) = 1.7925 \times 10^{- 3}$ for $ℏ ω_{o} = 2.6264$ eV for the data of Spicer, from which values of $R (ω)$ and $p (E)$ can be investigated.

FIG. 27.

View large Download slide

The measured photoemission data of Cs $_{3}$ Sb (Spicer and Taft data digitally extracted from Refs. 11 and 132, respectively) compared to Eq. (121) (red, “SMM”), and Eq. (122) (blue dashed, “Spicer”) for a bandgap $E_{g} = 1.6$ eV and electron affinity $E_{a} = 0.4$ eV.

3. Optical and scattering models

The laser penetration depth $δ (ω)$ and reflectivity $R (ω)$ are both dependent on the properties of the complex index of refraction $\hat{n} = n + i k$ for metals and semiconductors,^133–135 which can be theoretically described using a Lorentz + Drude + Resonant term (LDR) model,¹²⁶ or experimentally tabulated.^136–139 For present needs, it is assumed that tabulations of $n$ and $k$ are available across the IR to UV range 1.5 eV $< ℏ ω <$ 6.5 eV most common to photocathodes of interest to photoinjectors^140–142 for accelerators and free electron lasers.

In general, the total dielectric constant

K (ω) = \hat{ε} / ε_{0}

is the sum of free (Drude) and bound (Lorentz) components and can be separated into real

K_{r}

and imaginary

K_{i}

parts which are related by the Kramers–Kronig relations in descriptions of the optical properties of dielectrics.^143,144 A Drude–Zener model treats electrons subject to an oscillatory force and undergoing scattering events and performs well in the infrared and smaller frequencies for both metals and semiconductors but requires quantum mechanical modifications in the optical and higher frequencies.^134,145 The Lorentz–Drude (LD) model amends the Drude–Zener model by accounting for the presence of dampening forces¹⁴⁶ and thereby accounts for metallic-like behavior of simple metals like aluminum, but additional terms are needed for the description of metals and semiconductors used as photocathodes:⁸ resonances appear in the optical and UV spectrum where photocathodes are operated. The LDR model is then defined by¹⁴⁷

\begin{matrix} K (ω) = 1 + χ_{f} (ω) + χ_{b} (ω), \\ χ_{f} (ω) = - \frac{f_{d} ω_{p}^{2}}{ω (ω + i Γ_{d})}, \\ χ_{b} (ω) = - \sum_{j = 0}^{N_{b}} \frac{f_{j} ω_{p}^{2}}{ω^{2} - ω_{j}^{2} + i ω Γ_{j}}, \end{matrix}

(123)

where the Drude term is

χ_{f} (ω)

⁠, the Lorentz term is the

j = 0

part of

χ_{b} (ω)

⁠, and the resonant terms are the

j \in (1, 2, \dots, N_{b})

terms in the summation in

χ_{b}

⁠, where

N_{b}

can be made as large as needed for a desired accuracy. The terms

f_{j}

are dimensionless oscillator strengths, and

Γ_{j}

are damping rates. The term

ω_{p}

is the plasma frequency defined by

ω_{p}^{2} \equiv \frac{q^{2} ρ_{o}}{ε_{0} m} = \frac{16 π Q}{m} ρ_{o},

(124)

where the first form is in the usual SI notation, the latter in a form appropriate for the units of Sec. I A,

ρ_{o} = k_{F}^{3} / 3 π^{2}

is the number density for a zero temperature electron gas with a Fermi level of

μ = ℏ^{2} k_{F}^{2} / 2 m

⁠, and

k_{F}

is the Fermi momentum (the reader is cautioned not to confuse the “

k

” of the index of refraction

\hat{n} = n + i k

with the

k

denoting wave number as in

k_{F}

⁠). As an example, consider copper, for which

μ = 7

eV, which entails that

k_{F} = 13.5546

nm

^{- 1}

⁠, and so

\begin{aligned} ℏ ω_{p} & = {[\frac{ℏ^{2}}{m} 16 π Q ρ_{o}]}^{1 / 2} = {[\frac{32 μ k_{F} Q}{3 π}]}^{1 / 2} = 10.7691 eV . \end{aligned}

(125)

In terms of them, the index of refraction terms can be evaluated⁸ via

\begin{matrix} 2 n^{2} = \sqrt{K_{r}^{2} + K_{i}^{2}} + K_{r}, \\ 2 k^{2} = \sqrt{K_{r}^{2} + K_{i}^{2}} - K_{r}, \end{matrix}

(126)

where

K = K_{r} + i K_{i}

⁠, from which

R (ω)

and

δ (ω)

are deduced from

\begin{matrix} R (ω) = \frac{{(n - 1)}^{2} + k^{2}}{{(n + 1)}^{2} + k^{2}}, \\ δ (ω) = \frac{c}{2 k ω} . \end{matrix}

(127)

Resonance terms can be parametrically added^133,148 and methods to incorporate as many resonance terms as needed to model metals like copper and semiconductors like Cs $_{3}$ Sb complete the Lorentz–Drude–Resonance (LDR) model.¹²⁶ Numerous additional terms can arise given the complexity of $K (ω)$ ⁠, which can make $N_{b}$ about 11 for copper, 18 for gold, 36 for Cs $_{3}$ Sb, and 56 for perovskites.

For purposes of illustration, an ad hoc low- $N_{b}$ model is illustrated here for copper, which has enough complexity to reveal the advantages of the method. The LDR technique, however, is better suited to complex materials like perovskites and multi-alkali antimonides, to allow the parameterization of density functional theory (DFT) calculations of complex materials to assess their utility, or to see the consequences of variations in stoichiometry to photoemissive materials.^32,126 The Drude components (⁠ $f_{d}, Γ_{d}$ ⁠) are found by considering the small $ω$ behavior of $K (ω)$ and the Lorentz components (⁠ $f_{0}, Γ_{0}, ω_{0}$ ⁠) by the large $ω$ behavior. The resonance terms are then found by removing the LD part from the measured data and accounting for spikes in $K (ω)$ through the introduction of Lorentzians characterized by (⁠ $f_{j}, Γ_{j}, ω_{j}$ ⁠). The outcome of the analysis results in a reasonable comparison even for a small $N_{b}$ ⁠, as shown in Fig. 28(a), for which $N_{b} = 5$ ⁠. The parameters are found as follows, where, for convenience, $γ_{j} \equiv ℏ Γ_{j}$ is used so that the more convenient photon energy $ℏ ω$ can be used:

The Drude terms are found from the small- $ω$ behavior of $K_{i}$ to be $f_{d} = 0.7$ and $γ_{d} = 0.0888$ eV.
Examination of $K_{i} (ω)$ near $ℏ ω_{0} = 45$ eV give $f_{0} = 32$ and $γ_{0} = 300$ eV.
The peak locations were set to ${ℏ ω_{j}} = (1.8, 2.4, 5, 14.5, 25)$ in [eV]; crude estimates of ${γ_{j}} = (1.5, 1.0, 4.0, 5.0, 5.0)$ in eV characterized the widths and ${f_{j}} = (- 0.05, 0.08, 0.56, 0.63, 0.38)$ the peak heights (observe the sign of the first term).

Once these factors are found, the LDR model can be used without further ado, and so the evaluation of

R (ω)

and

δ (ω)

quickly accomplished thereafter for arbitrary

ℏ ω

⁠. Although the

γ_{j}

were ad hoc the

f_{j}

were set by hand as a pedagogical exercise here, a more rigorous search for the

γ_{j}

and an exact method to find

f_{j}

by inverting a matrix equation treating them as a vector to be solved for provides a much better correspondence.¹²⁶ It is seen that even a small

N_{b} = 5

number of

f_{j}

is sufficient to give reasonable

K_{t h}

lines for purposes of simulation. Still, the difference between the LD and LDR models is shown in Fig. 28(b), revealing the importance of the resonances for the visible and UV frequencies that photocathodes are excited by (Fig. 29).

FIG. 28.

View large Download slide

(a) Comparison of measured dielectric data for copper (symbols) to the theoretical LDR model with a small number $N_{b}$ of ad hoc values of $(f_{j}, γ_{j}, ℏ ω_{j}$ ⁠). (b) Comparison of the LD model to the LDR model including (⁠ ${j} \in (1, \dots, N_{b})$ ⁠) resonance terms.

FIG. 29.

$(a) Evaluation using Eq. (126) of index of refraction terms ( n , k ) for LDR parameters shown in Fig. 28. (b) Evaluation using Eq. (127) of reflectivity and penetration depth ( R , δ ) for LDR parameters shown in Fig. 28, where for visualization, δ ( ω ) / 20 is shown.$

View large Download slide

(a) Evaluation using Eq. (126) of index of refraction terms $(n, k)$ for LDR parameters shown in Fig. 28. (b) Evaluation using Eq. (127) of reflectivity and penetration depth $(R, δ)$ for LDR parameters shown in Fig. 28, where for visualization, $δ (ω) / 20$ is shown.