Ab Initio calculation of field emission from metal surfaces with atomic--scale defects

TL;DR

This paper combines density functional theory and quantum transport calculations to analyze how atomic-scale defects affect the work function and field emission from metal surfaces, revealing local enhancements due to defects.

Contribution

It introduces a general methodology for calculating field emission currents from nano-featured metal surfaces, specifically applied to defected Cu(111) surfaces.

Findings

Defects can locally increase emission current density.

The increase is mainly due to work function reduction, not field enhancement.

Fowler--Nordheim equation remains valid with adjusted work function.

Abstract

In this work we combine density functional theory and quantum transport calculations to study the influence of atomic--scale defects on the work function and field emission characteristics of metal surfaces. We develop a general methodology for the calculation of the field emitted current density from nano-featured surfaces, which is then used to study specific defects on a Cu(111) surface. Our results show that the inclusion of a defect can significantly locally enhance the field emitted current density. However, this increase is attributed solely to the decrease of the work function due to the defect, with the effective field enhancement being minute. Finally, the Fowler--Nordheim equation is found to be valid when the modified value for the work function is used, with only an approximately constant factor separating the computed currents from those predicted by the Fowler--Nordheim equation.

Tables

Table 1. Table 1: The work functions determined for the different surface defects with converged surface sizes. Δ ​ ϕ Δ italic-ϕ \mathinner{\mathrm{\Delta}\phi} is the difference to the clean surface work function.

System	Work function $\varphi \left(\mathrm{eV}\right)$	$\Delta \varphi \left(\mathrm{eV}\right)$	$\Delta \varphi \left(\%\right)$
Clean surface	$4.76$	$0.00$	$0.0$
Step defect	$4.66$	$-0.10$	$-2.1$
Adatom defect	$4.44$	$-0.32$	$-6.7$
Pyramid defect	$4.25$	$-0.51$	$-10.7$

Table 2. Table 2: Field emission electron currents for the different systems and electric fields at zero temperature in A   m − 2 times ampere meter 2 \mathrm{A}\text{\,}{\mathrm{m}}^{-2} . The results for the Schottky–Nordheim barrier are shown for comparison.

	$1\mathrm{GV}{m}^{-1}$		$2\mathrm{GV}{m}^{-1}$		$3\mathrm{GV}{m}^{-1}$
S.–N. barrier	$3.96$	$\times {10}^{-17}$	$2.19$	$\times {10}^{-1}$	$4.69$	$\times {10}^4$
Clean surface	$8.75$	$\times {10}^{-18}$	$9.23$	$\times {10}^{-2}$	$1.76$	$\times {10}^4$
Step defect	$1.33$	$\times {10}^{-16}$	$1.52$	$\times {10}^{-1}$	$1.84$	$\times {10}^4$
Adatom defect	$1.15$	$\times {10}^{-14}$	$1.21$		$7.98$	$\times {10}^4$
Pyramid defect	$9.53$	$\times {10}^{-12}$	$2.66$	$\times {10}^1$	$7.13$	$\times {10}^5$

Equations

$$J(\varphi ,F,T)=e{\int }_0^{\infty }\mathrm{d}{E}_{z}N(E_z,T)D(E_z,\varphi ,F),$$

$$J(\varphi ,F,T)=\frac{a{F}^2}{\varphi {\tau }^2(F)}\exp \left[-\frac{\nu (F)b{\varphi }^{3/2}}{F}\right],$$

$$\nu (F)=\frac{G(E_{\text{F}},\varphi ,F)}{G_{\text{ET}}(E_{\text{F}},\varphi ,F)}$$

$$J(\varphi ,F,T)=\lambda \frac{a{\beta }^2{F}^2}{\varphi {\tau }^2(\beta F)}\exp \left[-\frac{\nu (\beta F)b{\varphi }^{3/2}}{\beta F}\right].$$

$$N(E_z,T)={\int }_{E_z}^{\infty }\mathrm{d}E\frac{f_{\text{FD}}(E,T)\rho (E)}{\sqrt{8mE}}.$$

$$J=\frac{e}{\sqrt{2m}}{\int }_0^{\infty }\mathrm{d}E\sqrt{E}{f}_{\text{FD}}(E)\rho (E)⟨D(E,\theta ,\varphi )⟩,$$

$$⟨D(E,\theta ,\varphi )⟩=\frac{1}{2\pi }{\int }_0^{2\pi }\mathrm{d}\varphi {\int }_0^{\pi /2}\mathrm{d}\theta \sin \theta \cos \theta D(E,\theta ,\varphi )$$

$$V_{\text{xc}}^{\text{im}}(z)=\frac{e^2}{4\pi {ϵ}_0}\frac{1}{4(z-z_{\text{im}})},$$

$$V_{\text{xc}}=f(x)V_{\text{xc}}^{\text{PBE}}+[1-f(x)]V_{\text{xc}}^{\text{im}}$$

$$f(x)=\{\begin{array}{ll}1.0,& x\le 0\\ {\mathrm{e}}^{-x/{\lambda }_x}+(1-{\mathrm{e}}^{-x/{\lambda }_x}){\mathrm{e}}^{-x},& x>0\end{array}$$

$$\underset{Y}{\underbrace{\ln \left(\frac{J}{F^2}\frac{\varphi }{a}\right)}}=-\frac{1}{\beta }\underset{X}{\underbrace{\frac{{\nu }_{0}b{\varphi }^{3/2}}{F}}}-\underset{\text{constant}}{\underbrace{\frac{{\nu }_{1}b{\varphi }^{3/2}}{\beta }+\ln \left(\lambda \frac{{\beta }^2}{{\tau }_0^2}\right)}}.$$

$${\beta }_{\text{app}}={\left(\frac{{\varphi }_{\text{a}}}{\varphi }\right)}^{3/2}\beta ,$$

$$V(z)=A-eFz-\frac{e^2}{4\pi {ϵ}_0}\frac{1}{4(z-z_{\text{im}})},$$

$$J=\frac{e}{\sqrt{2m}}{\int }_0^{\infty }\mathrm{d}E\sqrt{E}{f}_{\text{FD}}(E)\rho (E)⟨D(E,\theta ,\varphi )⟩$$

$$⟨D(E,\theta ,\varphi )⟩=\frac{1}{2\pi }{\int }_0^{2\pi }\mathrm{d}\varphi {\int }_0^{\pi /2}\mathrm{d}\theta \sin \theta \cos \theta D(E,\theta ,\varphi )$$

Full text

I Introduction

Field electron emission (FE) from metal surfaces plays a crucial role in various aspects of modern technology. It is widely exploited in various devices that require cold electron sources De Heer et al. (1995); Saito and Uemura (2000); Okano et al. (1996); Jensen (1999), such as electron microscopes Crewe et al. (1970), flat displays Liu et al. (2012); Kuznetzov et al. (2010), electron beam lithography Vieu et al. (2000), electrospray propulsion Lemmer (2017), etc Egorov and Sheshin (2017). On the other hand, it often appears parasitically (in this case, it is usually referred to as dark current), limiting the performance or even igniting catastrophic vacuum breakdown events in a wide range of devices that involve metal surfaces exposed to high electric fields. Such devices span from micro- and nano-electronic capacitors Lyon and Hubler (2013); Ducharme (2009), macroscopic high-voltage devices such as vacuum interrupters and X-ray tubes Latham (1995); Slade (2007), up to large-scale apparatuses such as fusion reactors Jüttner (2001); McCracken (1980) and particle accelerators Dobert et al. (2005); Kim et al. (2008); Aicheler et al. (2012); Boland et al. (2016); Nanni et al. (2015).

Vacuum breakdowns (also known as vacuum arcs) are initiated at sites on the metal surface which exhibit strong field electron emission (Latham, 1995; Anders, 2008), even when the applied surface electric fields are in the range $30250$ $\text{\,}\mathrm{MV}\text{\,}{\mathrm{m}}^{-1}$, far lower than the FE regime ($110$ $\text{\,}\mathrm{GV}\text{\,}{\mathrm{m}}^{-1}$). The systematic experimental study of such emitting sites Beukema (1972); Descoeudres et al. (2009a, b) has shown an apparent field enhancement factor $\beta$ in the range $20140$. In other words, the current-voltage characteristics of such sites are consistent with a local field $\beta$ times higher than the applied macroscopic one. The most common interpretation of this apparent field enhancement is geometric, i.e. explained by assuming the existence of microscopic protrusions on the metal surface that locally enhance the electric field Beukema (1972); Descoeudres et al. (2009a, b). However, $\beta$ values of the aforementioned order correspond to protrusions with an aspect ratio in the same range, which have never been observed on the studied surfaces, even when high-resolution electron microscopy is utilized Kildemo et al. (2004).

The above discrepancy indicates that the classical Fowler–Nordheim (FN) theory Fowler and Nordheim (1928); Nordheim (1928); Murphy and Good (1956) with its standard interpretation for analyzing experimental data Forbes (2009) might be insufficient to describe FE under certain conditions. The FN equation Murphy and Good (1956) that canonically describes FE was derived assuming a mathematically flat planar metal surface with a given work function, a free–electron model for the metal bulk and a semi–classical approximation for the tunneling probability. It has already been shown that the FN equation is valid only for surfaces with radii of curvature larger than $\approx$20\text{,}\mathrm{nm}$$ Cutler et al. (1993); He et al. (1991). Although the behaviour of non-planar emitters has been studied extensively Zuber et al. (2002); Edgcombe and de Jonge (2007); Cabrera et al. (2013); Kyritsakis and Xanthakis (2015, 2016); Kyritsakis and Djurabekova (2017), these works relied on analytical approximations for the surface potential barrier, the assumption of a smooth, well defined surface for the emitter and a semi-classical approach for the transmission probability.

One mechanism that could possibly contribute in the explanation of the above unnaturally high apparent field enhancement factors is the local increase of the field emission current due to atomic–level surface defects such as adatoms, adsorbates, atomic steps etc. Such defects have already been shown to have a significant effect on the work function Djurabekova et al. (2013). Furthermore, field electron emission from atomic–scale features cannot be conceptually treated with the classical FN theory, since most of its assumptions are not valid. However, ab initio computational techniques Neugebauer and Scheffler (1992); Lepetit (2017); Groth et al. (2014) offer the ability to overcome the above problems and directly calculate the emission current densities for any surface configuration, without adopting most of the mathematical simplifications of the classical theory.

In this paper, we develop an ab initio computational method to calculate the field emitted current from defective metal surfaces. Our method combines density functional theory (DFT) electronic structure calculations with quantum transport methods. We use our method to obtain the field emission characteristics from defective Cu surfaces, with three different types of features, namely adatoms, mono-atomic steps and two-layer pyramids. We find that the current density and the corresponding apparent field enhancement increase on the defective surfaces, compared to the clean one. However, this increase is attributed mainly to the decrease of the corresponding work function. The methodology developed here can be used in the future to perform similar calculations for various kinds of surface defects such as adsorbates, oxide layers and larger nano-features.

The structure of the paper is the following: Section II will introduce the theoretical background for this work. Section III will describe the calculation of the field emission current. Section IV will show the results of the calculations and their interpretation using Fowler–Nordheim plots. Limitations of the methodology and the interpretation of the results will be discussed in section V. Finally, the conclusions drawn from the work will be given in section VI.

II Theory

II.1 Field emission theory

Field emission is still most often described using the Fowler–Nordheim equation developed in the 1920s Nordheim (1927); Fowler and Nordheim (1928); Nordheim (1928). The Fowler–Nordheim theory makes several approximations in order for the field emission current to be analytically solvable. The metal is assumed to be described by the free–electron jellium model, with the ionic potential having a discontinuity at the surface, where it changes to the vacuum potential. The surface is assumed to be perfectly flat, reducing the problem to one dimension. The jellium potential is smoothed by the image potential, which has its image plane at the jellium surface.

Murphy and Good Murphy and Good (1956) give a rigorous mathematical treatment of the problem with the assumptions above, giving the emission current as an integral

[TABLE]

where $E_{z}=\hbar^{2}k_{z}^{2}/(2m)$ is the normal (to the emitting surface) energy component, $\bm{\mathrm{k}}=(k_{x},k_{y},k_{z})$ is the wave vector, $m$ is the electron rest mass, $\phi$ is the work function, $F$ is the applied electric field on the surface and $T$ is the temperature. $N(E_{z},T)$ is the supply function, which gives the normal energy–dependent flux of electrons to the surface, and $D(E_{z},\phi,F)$ is the transmission probability, the probability for an electron to tunnel through the surface potential barrier into the vacuum. Using the above assumptions along with the semi-classical JWKB approximation Landau and Lifshitz (1977) for the transmission coefficient, the integral can be solved analytically for low temperatures. For $T=0$ it yields the Fowler–Nordheim equation

[TABLE]

where $a$ and $b$ are universal constants and $\nu(F)$ and $\tau(F)$ are well-known mathematical functions. Murphy and Good also give a correction factor to the Fowler–Nordheim equation for low non-zero temperatures.

Forbes Forbes (2004); Forbes and Deane (2007); *Forbes2007corr gives interpretations of $\nu$ and $\tau$, linking them to the Gamow exponent $G=-\ln D(E_{z},\phi,F)$ and its derivative with respect to $E_{z}$ at the Fermi level. The barrier enhancement factor

[TABLE]

is the ratio of the Gamow exponent at the Fermi level to that of the corresponding exact triangular barrier, while $\tau(F)$ is the same for the derivative of the Gamow exponent with respect to the normal energy.

Two additional effects are often added to eq. (2) to obtain the technically complete Fowler–Nordheim equation Fischer et al. (2013). Due to geometric effects not taken into account when assuming the planar surface, the electric field at the surface can be larger than the macroscopic applied electric field $F$. This is taken into account using the field enhancement factor $\beta$, replacing $F$ with $\beta F$. Also, a multiplicative constant $\lambda$ is added to the equation to include the effects of the band structure and other effects of real materials which do not appear in the jellium model. With these changes, the emission current density is given by Fischer et al. (2013)

[TABLE]

II.2 Calculation from first principles

II.2.1 Deviation from Fowler–Nordheim theory

In the current work, we abandon simplifying approximations used in the standard FN theory, in order to provide a more realistic description of the field emitting system from first principles.

First, the supply function $N(E_{z},T)$ is modified to allow for arbitrary density of states. In the standard FN theory which is based on the free electron model, the electron states are considered to be plane waves, with an isotropic parabolic dispersion relation $E=\hbar^{2}k^{2}/2m$ and a density of states that has a square root dependence on the total kinetic energy. In this work, we use the density of states as calculated by DFT, while approximating the Bloch states as plane waves with free electron dispersion. Thus, the supply function can be computed by integrating the density of occupied states over the slice of the reciprocal space which corresponds to a certain normal energy:

[TABLE]

In the above formula, $f_{\text{FD}}$ is the Fermi–Dirac occupation function, $\rho$ is the density of states and $E$ is the total electron kinetic energy. Additionally, although we do not demand zero temperature for the occupation function, the minute perturbative effect of the electron relaxation at non-zero temperatures on the potential is neglected.

Second, we step away from the simplified model of a jellium metal with a perfectly flat surface. Instead we model the metal surface in density functional theory as an atomic system for the computation of the potential barrier. However, the electron momentum components parallel to the surface were approximated as having no effect on the transmission probability, which is exact only under the assumption of a flat surface. This approximation is necessary in order to use equation (1). Otherwise, computing the differential current density for each electron kinetic energy would have required integrating the transmission probability over the isoenergy surface, which is too computationally demanding. 111The calculation replacing equation (1) would then take the form

$J=\frac{e}{\sqrt{2m}}\int_{0}^{\infty}\mathinner{\mathrm{d}E}\sqrt{E}f_{\text{FD}}(E)\rho(E)\left\langle D(E,\theta,\phi)\right\rangle\,,$

(6)

where

$\left\langle D(E,\theta,\phi)\right\rangle=\frac{1}{2\mathrm{\pi}}\int_{0}^{2\mathrm{\pi}}\mathinner{\mathrm{d}\phi}\int_{0}^{\mathrm{\pi}/2}\mathinner{\mathrm{d}\theta}\sin\theta\cos\theta D(E,\theta,\phi)$

(7)

is a weighted average of the transmission probability over the upper hemisphere of the isoenergy surface and $\theta$ and $\phi$ are the polar and azimuthal angles in $\bm{\mathrm{k}}$-space, respectively. Equation (1) with the supply function from equation (5) follows from this when assuming that $D(E,\theta,\phi)=D(E_{z})$, where $E_{z}=E\cos^{2}\theta$.

Finally, in order to calculate the transmission probabilities, we did not use the JWKB approximation Landau and Lifshitz (1977) as in the classical theory Murphy and Good (1956), or perturbative methods as in recent ab-initio approaches Lepetit (2017). We computed them numerically by using quantum transport theory, as described in section III.3.

II.2.2 Image potential

Since the potential obtained directly from DFT calculations can be trusted at the level allowed by the used functional, it is necessary to assess the reliability of the selected functional in all regions of interest, which in our case are bulk, surface and vacuum. Although the Perdew–Burke–Ernzerhof (PBE) functional Perdew et al. (1996); *Perdew1996err used in this work is accurate in and near the metal, it fails to reproduce the long-range exchange-correlation potential in the vacuum Eguiluz et al. (1992). The true exchange-correlation potential in the vacuum over a smooth metal surface is known to converge to the image potential Bardeen (1940); Harbola and Sahni (1993); Lang and Kohn (1973)

[TABLE]

where $z_{\text{im}}$ is the location of the electrical surface, also known as the image plane. The exchange-correlation potential computed using the PBE functional decays exponentially into the vacuum, which makes the computed potential barriers too high.

In order to improve the description of the potential barrier, the image potential was externally introduced as an approximation of the true exchange-correlation potential in the vacuum region. The image plane was assumed as being flat and parallel to the $x$-$y$ plane. Its location is known to be the $z$ coordinate of the center of mass of the charge layer induced by the applied electric field Chulkov et al. (1997, 1999); Serena et al. (1986).

Having computed the exchange-correlation potential $V_{\text{xc}}^{\text{PBE}}$ in the vacuum region with DFT, the total potential was modified as follows to take into account the image potential. A smooth transition

[TABLE]

from the exchange-correlation potential computed using the PBE functional to the image potential was introduced, where

[TABLE]

is the a weighting function and $x=(d-d_{0})/\lambda_{d}$ is a rescaling of the distance $d=\max(r,z^{\prime})-z^{\prime}_{\text{im}}$. $r$ is the distance of the point for which the potential is being determined from the nearest atom, $z^{\prime}$ is the distance of the point from the image plane and $z^{\prime}_{\text{im}}$ is the distance of the image plane from the topmost layer of atoms in the metal, not counting the defect atoms. The form of the smoothed Heaviside weighting function $f(x)$ ensures that the resulting exchange-correlation potential is continuous and smooth. The $d_{0}$ parameter determines how far from the image plane the transition begins, $\lambda_{d}$ determines how long the transition distance is and $\lambda_{x}$ determines how abrupt the beginning of the transition is.

The computation of the center of mass was done by splitting the metal slab into two parts in the $z$ direction so that the total charge on each side was that predicted by Gauss’ law for the surface charge, and then calculating the quotient of the excess dipole moment of each side and the corresponding surface charge. The excess dipole moment was determined by subtracting the dipole moment of the same surface with no applied electric field from the dipole moment of the surface being studied.

III Methods

III.1 Simulation setup

For this study we selected to to simulate the (111) orientation of a Cu surface, since the expected effect on the work function from atomic defects on this surface is the strongest compared to other orientations. Three different types of surface defects were studied along with the clean surface, all shown in figure 1.

The first defect type was a step defect consisting of an additional layer of FCC atoms covering half of the surface. To exclude any finite size effect of the step due to periodic boundaries in lateral directions, the length of the cell perpendicular to the step was varied. The second defect type was a single adatom in an FCC lattice site on the surface. With periodic boundary conditions this resulted in a surface with an adatom density given by the inverse of the area of the supercell. The surface was kept approximately square and its size was varied. The third defect was a two layer triangular pyramid with the lower layer atoms in FCC lattice sites. The system otherwise resembled the adatom defect system. The surface was kept approximately square and its size was varied.

We carried out all DFT calculations using the version 5.3.5 of the Vienna Ab initio Simulation Package (VASP) Kresse and Hafner (1993, 1994); Kresse and Furthmüller (1996a, b) with the Projector Augmented-Wave method (PAW) Blöchl (1994); Kresse and Joubert (1999). The used pseudopotentials were from the version 54 dataset supplied with VASP. The used pseudopotential for copper treated the eleven 3d and 4s electrons as valence electrons and the rest of the electrons as a frozen core.

III.1.1 Calculation details

All calculations used the Perdew–Burke–Ernzerhof (PBE) exchange-correlation functional Perdew et al. (1996); *Perdew1996err. The electronic convergence criteria was set to a change of less than ${10}^{-7}\text{\,}\mathrm{eV}$ in the total energy and the Kohn–Sham eigenvalues. The convergence criterion for the ionic relaxation was set to ${10}^{-6}\text{\,}\mathrm{eV}$ in the total energy, which resulted in remaining forces of order ${10}^{-3}\text{\,}\mathrm{eV}\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-1}$.

The convergence of the plane wave energy cutoff, $\bm{\mathrm{k}}$-point sampling, smearing method and parameter and lattice constant was tested using a bulk FCC copper system with one atom in the simulation cell. The chosen parameters were $E_{\text{cut}}=$800\text{,}\mathrm{eV} for the energy cutoff and $ka\geq$45\text{\,}\mathrm{\SIUnitSymbolAngstrom} for the $\bm{\mathrm{k}}$-point sampling, where $a$ is the simulation cell size and $k$ the number of $\bm{\mathrm{k}}$-samples in each direction. A value $\sigma=$0.2\text{,}\mathrm{eV} was chosen for the first order Methfessel–Paxton Methfessel and Paxton ([1989](#bib.bib65)) smearing parameter, and $a_{0}=$3.635\text{\,}\mathrm{\SIUnitSymbolAngstrom} for the lattice constant. The energy cutoff and $\bm{\mathrm{k}}$-point sampling were chosen so that the estimated error from each was less than $1\text{\,}\mathrm{meV}$, the smearing parameter was chosen to give a small energy difference to using the tetrahedron smearing method with Blöchl’s corrections Blöchl et al. (1994), and the lattice constant was determined with a resolution of $1\text{\,}\mathrm{m\SIUnitSymbolAngstrom}$.

The unit cell for the surface calculations was chosen to be the minimal orthorhombic FCC unit cell with the $(111)$ direction along the $z$ axis and the $x$ direction being along the shorter cell edge. The unit cell contains three layers of two atoms. The lowest layer has the atoms in the corner and in the middle of the cell, and the atoms in the other layers are shifted by a third of the cell length in the $y$ direction each. The ionic positions were relaxed in all surface calculations with the bottom $2$ layers being fixed to the bulk positions.

We also performed convergence tests with regard to the number of slab layers and the thickness of the vacuum layer, using a system consisting of a single unit cell in the $x$-$y$ plane with a variable number of slab layers and variable vacuum thickness. Additionally, we studied the convergence of the work function and the interlayer spacing of the slab. In these tests, we found that the inter-layer spacing converges for a slab thickness of $8$ atomic layers. A vacuum layer of approximately 30 Å was sufficient to achieve convergence of the work function within 10 meV uncertainty while leaving enough vacuum space for the defect structures.

The convergence with regard to the size of the supercell parallel to the surface was done separately for each defect type. The used criterion for convergence was the effect of the defect on the work function being approximately inversely proportional to the size of the surface. The criterion effectively means that the area furthest from the defect is equivalent to the clean surface. The resulting surface sizes were $1\times 8$ unit cells for the step defect and $4\times 2$ unit cells for the adatom and pyramid defects.

III.2 Potential barriers

The potential barrier of interest for the tunneling of electrons into the vacuum is the one found at the surface of the slab. The potential oscillations inside the slab due to the atoms are far below the Fermi level and not of interest in the current work. The total local potential was taken from the DFT calculations described above merged with the image potential as described in section II.2.2.

The parameters for the transition from the PBE potential to the image potential were found by optimizing them for the clean surface, ensuring that the potentials merge smoothly near the region where the distance between them is the shortest (see Fig. 2). This optimization resulted in choosing the following values for the parameters used in expression (10): $d_{0}=$1.5\text{,}\mathrm{\SIUnitSymbolAngstrom}, $\lambda_{d}=$1.0\text{\,}\mathrm{\SIUnitSymbolAngstrom} and $\lambda_{x}=$0.75\text{,}\mathrm{\SIUnitSymbolAngstrom}$$.

The center of mass of the excess charge was found to be $1.56\pm 0.08\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ above the top atomic plane in the slab with the clean surface. This is slightly further out than the result for jellium slabs of $\approx$1.13\text{,}\mathrm{\SIUnitSymbolAngstrom}$$ Chulkov et al. (1999).

The exchange-correlation potential $V_{\text{xc}}^{\text{PBE}}$ was computed as the difference between the total local potential and the Hartree potential, which are outputs of VASP. Due to unphysical noise in the exchange-correlation potential approximately $5\text{\,}\mathrm{\SIUnitSymbolAngstrom}8\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ above the surface layer, $V_{\text{xc}}^{\text{PBE}}$ was numerically smoothed. The smoothing removed the noise while only negligibly effecting the form of the potential. We performed the smoothing by applying to each grid line in the $z$ direction a seventh order, $71$-point ($\approx$1.75\text{,}\mathrm{\SIUnitSymbolAngstrom}$$) Savitzky–Golay filter. After that, we applied two uniform filters, the first with a $5\times 5\times 5$ grid point stencil and the second with a $3\times 3\times 3$ grid point stencil.

Figure 2 demonstrates the shape of the modified exchange-correlation potential in a solid line. In the same figure, the original PBE and image potentials are shown in dashed lines with green and red colors, respectively.

III.3 Quantum transport

This section describes the steps necessary to calculate the field emission current densities from the studied surfaces of interest (see Fig. 1), using the supply function and the potential grid from the DFT calculations, corrected to include the image potential.

The numerical quantum transport calculations were carried out using the Kwant code Groth et al. (2014) (version 1.3.2), which implements the computation of transport properties for systems described by tight-binding Hamiltonians. The latter result from applying the finite difference method on the single-particle Schrödinger equation, assuming a certain potential distribution $V$, which is pre-calculated by DFT. The resulting tight-binding Hamiltonian includes only onsite terms $V+2\sum_{i}T_{i}$ and nearest-neighbor interactions $-T_{i}$, where $T_{i}=\hbar^{2}/(2ma_{i}^{2})$ is the hopping energy for the grid spacing $a_{i}$ in the direction $i$. In the tight-binding formalism, semi-infinite leads with uniform constant potential are defined for the metal and vacuum regions perpendicular to the surface. Then the scattering matrix of the transport modes in the leads is calculated by the wave function matching method.

Since the quantum transport calculations were done in the jellium approximation, it was necessary to define the constant potential level inside the metal. The output from DFT only gives the potential barrier in relation to the Fermi level, so the constant potential in the metallic lead had to be determined otherwise. Here we adopted the standard practice of choosing it as the bottom of the valence band Forbes (2009); Jensen (2007). However, we found that the if the zero potential is sufficiently deep (several eV below $E_{F}$), its position does not affect significantly the resulting transmission coefficients. The density of states in all calculations was nonzero only in a continuous area around the Fermi level. Hence, our choice is determined by the lowest energy valence electron state, which has zero kinetic energy in the used electron model.

After adding the image potential in the vacuum region, the total potential had to be processed in order to be used in the quantum transport calculations. Details of this procedure are given in appendix A. The shape of the processed potential is shown in figure 3.

We used the processed potential as the scattering section of the system in the quantum transport calculations. Identical leads were attached to the metal and the vacuum. The leads were at the zero potential and had the same translational symmetry as the potential grid, so the wave functions in the leads were plane waves. While this is correct inside the metal in the jellium model, the potential far away in the vacuum should decay linearly due to the electric field. However, this difference in the potential is far below the Fermi level, so its effect on the calculations was negligible.

Using the planar surface approximation for the supply function means that the momenta of the incoming waves parallel to the surface are approximated as irrelevant for the transmission probability. They were therefore set to zero, making the only relevant incoming mode the one with zero transverse energy. The transmission probability was computed as the sum of the transmission probabilites from the lowest mode in the metal lead to all modes in the vacuum lead. The single-mode transmission probabilities were computed from the scattering matrix, the primary output of the Kwant calculation.

The transmission probabilities were computed for energies from $3\text{\,}\mathrm{eV}$ below the Fermi level to $5\text{\,}\mathrm{eV}$ above the Fermi level with a resolution of $0.1\text{\,}\mathrm{eV}$. The resolution was changed to $0.025\text{\,}\mathrm{eV}$ for $0.1\text{\,}\mathrm{eV}$ to either side of the Fermi level due to the importance of this area for the field emission current.

After computing the transmission probabilities, the differential current densities were computed as the product of the transmission probability and the supply function. The total current density from the studied surface was then computed by integrating the differential current density over the normal energy, as in eq. (1). The numerical integration was done using Simpson’s rule.

IV Results

IV.1 Work functions

We computed the work function for each surface as the difference between the converged Hartree potential in the vacuum and the Fermi level using the calculations without applied electric field. The Hartree potential was fully converged to the vacuum value in each calculation, as can be seen in figure 4.

The work function values for the converged surface sizes are shown in table 1. One can see that, as expected, the work functions are lower for the defect systems than for the clean surface. The effect on the work function is largest for the pyramid defect and smallest for the step defect, meaning that the defects which cause a rougher surface have smaller work functions. This is expected due to the effect of Smoluchowski smoothing Smoluchowski (1941).

The work function of the clean surface is expected to be lower than the experimental value of $4.94\text{\,}\mathrm{eV}$ Gartland et al. (1972), as the PBE functional underestimates work functions by $0.3\text{\,}\mathrm{eV}$ on average De Waele et al. (2016). The computed value is in good agreement with the literature value for the PBE functional, $4.78\text{\,}\mathrm{eV}$ Da Silva et al. (2006). The reduction of the work function by $6.7\text{\,}\mathrm{\char 37\relax}$ for the adatom defect is in good agreement with the calculations of Djurabekova et al. Djurabekova et al. (2013), which yielded a reduction of $5.9\text{\,}\mathrm{\char 37\relax}$ for the copper $(100)$ surface using the same criteria for the size of the surface.

The validity of comparing the computed work functions of different systems is discussed in section V.1.

IV.2 Transmission probabilities and emission currents

Figure 5 shows the computed transmission probabilities as a function of the normal energy. One can see that the transmission probability is higher for a system with a smaller work function or a larger applied electric field. The order of the transmission probabilities for the different defect types is the same for all applied electric fields, and is the same as in table 1.

All transmission probability curves are qualitatively similar to those for the Schottky–Nordheim barrier. Notably, all of the curves are almost parallel near the Fermi level, with differences in their slope appearing only at higher energies near the top of the barrier, which are irrelevant for field emission due to the Fermi–Dirac statistics causing the supply function to vanish. The curves are steeper for the lower applied fields, because the width of the barrier increases faster as the energy decreases.

The curvature of the plots near the Fermi level is small, which justifies the use of the standard Fowler–Nordheim equation for finding the field enhancement factors, as its derivation is based on a linear approximation of the transmission probability near the Fermi level.

The computed densities of states for the different systems differed from each other only negligibly. An example can be seen in figure 6. They closely resembled the density of states of bulk copper, with a small background from the 4s electrons and a peak approximately $1.5\text{\,}\mathrm{eV}$ below the Fermi level from the 3d electrons. Since the energies with large densities of states are significantly below the Fermi level, they are not of significant importance for field emission. Neglecting the peak, one could approximate the density of states by the free electron parabola, which would result in the introduction of a constant correction factor to the emission current.

The differential current densities $j(E_{z})$, the emission current density per unit normal energy from electrons with the normal energy $E_{z}$, is shown at room temperature as a function of the normal energy in figure 7. The temperature determines the slope of the decrease above the Fermi level, which becomes infinite for zero temperature. Thus, the finite temperature contributes to the emission current from energies above the Fermi level.

The qualitative comparison of the curves for different systems similar to the one for the transmission probabilities, since the density of states which also affects $j(E_{z})$ is almost the same for all systems. However, the irregular shape of the computed density of states has a visible effect on the $j(E_{z})$ curves, appearing as a slight increase of $j(E_{z})$ at $E_{z}\approx$ $-1.5\text{\,}\mathrm{eV}$.

The integrated total currents are listed in table 2. One can see that the current density for the Schottky–Nordheim barrier is larger than the current density computed for the clean surface by a factor of approximately $25$ for all computed electric field strengths. Overall, one can see that the differences between the systems become smaller as the electric field increases. The increase of the emission current caused by the defects is large, about $3$ orders of magnitude for the single adatom and $6$ for the pyramid defect with the smallest electric field.

IV.3 Fowler–Nordheim plot analysis

In order to assess the effects of the various surface irregularities on the parameters usually extracted by field emission measurements, we shall perform Fowler–Nordheim plot analysis on the calculated current densities. For this purpose, we shall linearize eq. (4) by approximating the parameter $\nu(F)$ as a linear function of the applied electric field, $\nu=\nu_{0}+\nu_{1}F$ and the parameter $\tau$ as a constant $\tau_{0}$. The parameters $\nu_{0}$ and $\nu_{1}$ were computed from least squares fits to the computed emission currents, and it was found that $\nu_{0}\approx 1$ for all systems while $\nu_{1}$ varied. The linear approximation was justified post hoc, as it leads to linear Fowler–Nordheim plots, which were observed. The resulting linearized form of the Fowler–Nordheim plots is

[TABLE]

The field emission currents computed for the different defect type and applied electric fields are plotted in a Fowler–Nordheim type plot of the form (11) in figure 8. One can see that the plots are linear, with slopes very similar to the F-N equation (S-N barrier), which is unity. We shall now define the effective field enhancement factor (or slope correction factor) $\beta$ for each system, as the inverse of the slope of each curve. We call it effective field enhancement because it accounts for all the possible effects that induce a change in the slope of the F-N plot, equivalent to an actual electrostatic field enhancement of the same magnitude.

In figure 8, curves have similar slopes that deviate less than 2% from unity. Thus, the effective field enhancement factors induced by the simulated defects are insignificant. This was also verified using linear fits to the plots, which show that $\beta$ is unity to the accuracy allowed by this work. This means that the increase in the emitted currents due to the surface defects can only be a result of the decrease of the local work function.

The remaining difference between the computed plots and the reference plot for the S–N barrier is simply an approximately constant factor, which is smaller for the studied systems than for the reference system. This means that the technically complete Fowler–Nordheim equation (4) is sufficient to describe the field emission from the studied surfaces, given that correct values for the work function $\phi$ and the constant pre-factor $\lambda$ have been obtained.

However, in standard experimental circumstances, it is impractical to obtain an accurate value of $\phi$ for defective surfaces, and a conventional clean-surface value is usually assumed. In this case, the change in the slope of the F-N plot due to the different $\phi$ is interpreted as field enhancement. The apparent field enhancement factor $\beta_{\text{app}}$, i.e. the $\beta$ extracted from the F-N plot when a certain value $\phi_{\text{a}}$ is assumed for the work function, can be deduced from equation (11) (neglecting small differences in $\nu_{0}$) as

[TABLE]

where $\phi$ is the real work function and and $\beta$ the effective field enhancement factor of the surface (unity in our calculations). The ratio of the assumed and real work functions is greater than unity because of the surface defects.

Figure 9 shows the same plots as figure 8, with the difference that the work function and barrier enhancement parameter are assumed to be those of the clean surface for all surface defects. One can see that the slopes of the plots are now different due to the scaling of the abscissa with the work function. The apparent field enhancement factor of the adatom and pyramid defects, as extracted using eq. (11) are now approximately $1.14$ and $1.24$ respectively.

V Discussion

V.1 Work function determination

Much of the analysis in this work is based on quantifying the effect of a surface defect on the work function of a metal surface. The work function as computed here, is effectively a property averaged over the whole cell surface. It is therefore clear that the effect of a single defect on the calculated work function depends on the size of the cell, or in other words the coverage density of defects. This means that the comparison of the work function changes among different defect types is meaningful only if they all have the same size, i.e. the same surface coverage.

Furthermore, the cell size has to be sufficiently large, so that the defects do not interact with their periodic images. In order to assure this, we adopted a method given in an earlier work by Djurabekova et al Djurabekova et al. (2013), which studied the effect of an adatoms and step defects on a copper $(100)$ surface on the work function. The change of the work function resulting from a change in the lateral size of the slab is studied and convergence is reached when the change in the work function due to the defect is inversely proportional to the surface size, $\mathinner{\mathrm{\Delta}\phi}\propto 1/S$, or equivalently proportional to the coverage. This corresponds to a situation where the local changes around the defect do not affect the area near the periodic boundaries and the defects can therefore be considered as isolated. For one dimensional defects such as the step, the convergence criterion must be modified to take into account that the coverage depends only on the length of the supercell perpendicular to the step, $\mathinner{\mathrm{\Delta}\phi}\propto 1/L_{y}$.

V.2 Apparent field enhancement

One important result of this paper is that small surface defects, such as the ones simulated here, do not cause any significant change to the slope of the Fowler–Nordheim plot, apart from the one attributed to the modification of the work function. In other words, they do not induce any significant effective field enhancement in the field emission characteristics of a metal surface. Therefore, such defects alone, cannot explain the increased field emission currents and apparent field enhancements measured experimentally from clean Cu surfaces Kildemo et al. (2004). However, the reduction of the work function, which is not taken into account when analyzing F–N plots (a uniform value $\phi=$4.5 eV is usually considered for Cu Kildemo et al. (2004)), might contribute in the overestimation of the field enhancement.

Although these apparent field enhancement factors due to work function decrease are small, they are produced by $\mathrm{\SIUnitSymbolAngstrom}$-scale defects of copper. Surface contaminants which lower the work function (see e.g. ref. Leung et al., 2003) and larger defects which cause more surface roughness can cause much larger apparent field enhancement factors, which might provide an explanation for at least part of the field enhancement factors reported in experimental works.

The effect of finding large field enhancement factors due to disregarding the local work function has been studied before Chen et al. (2012); Qian et al. (2012), where it was found that the field enhancement factors for clean surfaces are approximately $15$ if the local work function is significantly lowered.

V.3 Limitations of the method

The method developed in this work is valuable for calculating field emission characteristics beyond the classical models. However, there are certain limitations imposed by the approximations that have been adopted.

First, the electrical surface of the slab, used to define the image plane, was assumed to be flat for all of the studied systems, even those with defects. This is due to the image potential only being known as the correct asymptotic exchange-correlation potential in the vacuum for a flat surface. This approximation limits the sizes of the surface and the defect structure which can be studied using the method described. A very high defect such as a nanorod would raise the computed electrical surface high above the surface of the surrounding clean surface. The effect of the image potential on the potential barrier there would therefore be removed and the potential would be too high. Also, using the image potential in the area around the nanorod would be questionable.

Second, the electric fields accessible by the specific DFT method used here are limit to about 3 GV/m. The standard field range for which field emission measurement from materials with Cu-like work function are taken is 3-7 GV/m Cabrera et al. (2013). Thus, the highest field achieved in our calculations lies in the lower limit for experimentally detectable field emission currents. However, although the current densities calculated for lower fields (1-2 GV/m) are not experimentally measurable, they are theoretically invaluable for calculating the slope of F–N plots and comparing with classical emission models.

VI Conclusions

We developed a methodology for performing ab initio calculations of field emission currents from metal surfaces using density functional theory and quantum transport calculations. The method was then used to analyze the effect of atomic-scale defects on copper surfaces on the field emission current and to determine the validity of the Fowler–Nordheim equation for these surfaces. The defects were found to decrease the work function, with the decrease being larger for defects with higher roughness. The largest local decrease of the work function was found to be $0.51\text{\,}\mathrm{eV}$, corresponding to an apparent field enhancement of 1.24 if the flat-surface work function were used. The analysis of the field emission currents determined using the developed method showed that the Fowler–Nordheim equation provides a good qualitative description of the field emission from surfaces with atomic-scale defects. For a quantitative comparison, the equation must be extended to the technically complete Fowler–Nordheim equation with a multiplicative pre-factor which depends on the specific material and surface configuration and can be calculated using our method. Building up on this work, the developed methodology can be used and expanded to calculate field emission characteristics of more complicated small surface structures and defects, or materials that cannot be described with the classical emission models.

Acknowledgments

H. Toijala and A. Kyritsakis were supported by the CERN K-contract (No. 47207461). The work of K. Eimre and V. Zadin was supported by the Estonian Research Council grant PUT 1372. F. Djurabekova acknowledges gratefully the financial support of the Academy of Finland (Grant No. 269696). We also acknowledge the grants of computer capacity from the Finnish Grid and Cloud Infrastructure (persistent identifier urn:nbn:fi:research-infras-2016072533).

Appendix A Details of processing the potential

After adding the image potential in the vacuum region, the potential still had to be processed in order to be used in the quantum transport calculations where the jellium model for the slab is considered. First, the potential was cropped to the dipole correction, removing the periodic image of the lower side of the slab. Second, most of the slab was cropped from the potential, leaving only $4\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ of slab below the top atomic plane of the substrate. This was only done to reduce the computational expense of the calculations and has no effect on the results.

Third, the potential in the vacuum was extended to the zero potential level of the jellium model. This was necessary because of the limited vacuum height used in the DFT calculations. A least squares fit was made to the potential for each grid line in the $z$ direction within the upper third of the vacuum region. The fitted function was of the form

[TABLE]

where $A$ is a fitting parameter, $F$ is the applied electric field and $z_{\text{im}}$ is the location of the image plane. Using the applied electric field instead of making it a fitting parameter is justified, as checking the potential output by VASP revealed that the electric field was equal to the applied electric field everywhere in this region within the margin of numerical error. Fitting the electric field would have caused slightly different values for neighboring grid lines and therefore discontinuities in the continued potential far from the surface.

Fourth, the potential grid was undersampled $12$ times, in order to reduce the required computational resources to feasible levels, while maintaining a sufficient accuracy. This was found to increase the computed transmission probability near the Fermi level by approximately $20\text{\,}\mathrm{\char 37\relax}$ compared to undersampling the potential only $4$ times (computing a reference using no undersampling was not feasible). The undersampling did not break the translational symmetry of the potential grid. The resulting potential grid had a spacing of approximately $0.32\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, allowing plane waves with up to $37.6\text{\,}\mathrm{eV}$ kinetic energy to be described, while the largest energy required for the calculations was approximately $15\text{\,}\mathrm{eV}$.

Finally, the potential inside the slab was clamped to the zero potential in order to emulate a jellium slab. For each grid line in the $z$ direction, the last place where the potential crossed the zero potential before the surface barrier was found. If this did not happen for a grid line, the $z$ coordinate of the top atomic plane was chosen. The potential in the area below the surface determined above was set to the zero potential.