A methodology of modeling nonplanar surfaces, in which the microscale features of the emission sites can be orders of magnitude smaller than the mesoscale features defining the active emission area, has been developed and applied to both ordered arrays of identical emitters and random variations characteristic of a roughened surface. The methodology combines a general thermal-field-photoemission model for electron emission, a point charge model for the evaluation of field enhancement factors and surface geometry, and a Ballistic-Impulse model to account for the trajectories of electrons close to the cathode surface. How microscale and mesoscale features can both undermine the estimation of thermal-field emission parameters, such as characteristic field enhancement and total current predictions, as well as give rise to changes in the distribution of transverse velocity components used to estimate beam quality features such as emittance that are important to photocathodes, is quantified. The methodology is designed to enable both the proper characterization of emitters based on experimental current-voltage data and the development of a unit cell model of emission regions that will ease the emission model demands in beam optics codes.

Mesoscale and microscale features coexist on thermal, field, and photoemission sources, either by design (as for arrays of emitters) or as a result of material properties (as for physical surface roughness or crystal face work function variation). They complicate material parameter characterization, determination of beam quality parameters, and predictions of cathode and device performance. The complications arise from how surface structure in the form of geometric shielding effects or emission patchiness affects thermal-field and photoemission properties and electron trajectories at three typical length scales: (i) microscale for feature sizes from 0.5 to 50 nm, (ii) mesoscale for feature sizes from 0.05 to 10 μm, and (iii) macroscale for features of a generally visible size. Previously, the relations between them were governed by the size of the emitter and its apex radius,1 but for feature sizes considered below, the dimensioned relations are approximately the same. The prediction of emission properties,2 the inference of microscopic structure and surface properties from measured current from structured materials,3 or the onset of breakdown4 is hampered by variations in both geometry (“geometric” surface roughness or “field enhancement”) and work function (“chemical” surface roughness, alternately called “patch”5–7 or “island effect”8 for polycrystalline surfaces or those with variations in coatings or hydrogen termination). Except for idealized configurations,9–14 exploration of the transition from field-dominated to thermal-dominated effects is difficult for a single emitter and generally not possible for an array of them using numerical methods that grid the surface. Photoemission from cathodes with surface roughness encounters analogous problems with respect to thermal and field effects on quantum yield15–19 and mean transverse energy (MTE) predictions.10,20,21

The canonical equations governing thermal, field, and photoemission [the Richardson-Laue-Dushman (RLD) equation, the Fowler-Nordheim (FN) equation, and the Fowler-DuBridge (FD) equation, respectively] can be obtained as limits of a general equation22–24 and, therefore, have ranges of applicability. Although the RLD and FN equations, in the form of their general usage,25,26 appear visually similar, a seemingly reasonable hypothetical “thermal-field” equation proposed by Millikan and Lauritsen27 was in fact errant.24 More forebodingly, warm field emission and tepid thermal emission conditions can coexist on the same emitter structure, so that properly predicting the current over a multidimensional emitter or extracting characteristic parameters is error-prone if using only the canonical equations. The difficulty is epitomized by monitoring the energy distribution as conditions evolve from thermal to field emission (e.g., Fig. 2 of Ref. 28, where the peak of the distribution shifts from the barrier maximum to near the Fermi level as the field increases).

Such behavior affects the inference of work function Φ, field enhancement β, and notional emission area from current-voltage or I(V) data and further has an impact on predictions of a beam’s MTE. This collection of characteristics will be referred to as the “emission parameters” below. The Point Charge Model (PCM),23,24 which represents an emitting protrusion by a linear array of point charges, enables quantifying them because the potential everywhere is analytic, and so surface fields can be accurately found even for microscale features on mesoscale structures. Analogous methods, such as the Line Charge Model (LCM)1,29 and the ring charge model,30 rely on similar strategies, although only the PCM and LCM are extensible to arrays or random placement of bumps as in surface roughness as will be done below. That flexibility is because (to paraphrase Eyges31) charges and only charges produce a potential, and if they produce an equipotential surface that matches an emitter surface, the solution is unique.

In the first part of the present study, the impact of array geometry for conical emitters is considered by characterizing how the pitch (or tip to tip separation) affects the field enhancement (β) and apex radius (an) factors. Using the PCM to construct a hyperbolic or ellipsoidal parameterization enables utilizing findings from prolate spheroidal models11,24,32–35 even for arrays. The transition from field-dominated emission to thermal dominated emission along the same emitter surface uses a general thermal-field equation for current density, allowing the total current from the emitter to be found. How typically unobservable microscale effects affect I(V) and the estimation of the emission parameters in the PCM is related to the multiplicative nature of field enhancement (aka “Schottky’s conjecture”13,36,37), and the predictive accuracy of conventional methods to estimate the emission parameters can be assessed when thermal and field mechanisms operate simultaneously.

In the second part of the present study, a randomized PCM is used to approximate surface roughness. The evaluation of local field enhancement and the specification of the surface elements become more difficult compared to the more orderly arrays of the first part and complicates how trajectories from the surface are found. The consequences of surface roughness on yield and MTE are of particular concern for photocathodes15,16,19,38–43 and serves as the impetus to find trajectories from the surface,21 but insofar as the FD equation allows the estimation of work function,44 the methods of the first section pertain in assessing thermal-field contributions. An analytic ballistic-impulse model21 is improved and applied to treating randomized surface roughness and allows roughness contributions to emittance for photocathodes to be evaluated.

This section develops techniques to estimate current-voltage characteristics from an ordered array of field-enhancing emitter tips with sufficiently low work function to also show appreciable thermal-field and pure thermal emission. This treatment is motivated by the possibility of coating the relatively new low work function material 12CaO7Al2O3:4e, or C12A7 electride,45 onto a field-enhancing patterned substrate for increased thermal-field emission. C12A7 is a binary eutectic in the calcium aluminate system notable for a unit cell consisting of several positively charged nanometer-scale “cages.” The unit cell balances charge by incorporating extra-framework anions such as O2, but after oxygen reduction, it supports free electrons and becomes a metallic conductor with work function 2.4 eV.46 The material has demonstrated A/cm2-level emission at several amperes and temperatures <1000 °C in hollow cathodes used for spacecraft plasma propulsion,47 albeit with limited lifetime, and is more commonly used commercially as an electron injection layer and a phosphor layer in organic light-emitting diodes (OLEDs).48 

The idea of coating a low work function or negative electron affinity material onto a field-emitting surface is not new; for example, diamond-coated silicon spikes and rare earth oxide-coated carbon nanotubes have been reported.49,50 However, C12A7’s combined low work function and nanoscale cage features have previously shown strong thermal-field emission even on relatively smooth (e.g., not explicitly patterned) samples.51,52 This section’s calculations will show that a low work function coating of C12A7 on micrometer-scale conical spike arrays with a modest apex radius of 50nm should substantially increase the field and thermal-field contributions to the total emission. These arrays are modeled on widely commercially available patterned sapphire substrates used in GaN LEDs to facilitate future experimental validation. Additionally, over wide ranges of input parameters, neither FN nor RLD treatments accurately capture the total emission from the coated array, falling orders of magnitude short in some cases.

The array model is based on a physical picture of a thin film C12A7 electride coating on a patterned sapphire substrate (PSS), as shown in Fig. 1. The individual sapphire tips are approximately conical with a height of 1.6μm, a base diameter of 2.60μm, and an apex radius of 55 nm. The sides of the cone are not linear but rather rounded with a radius of curvature approximately 3× the emitter height. At present, a multiplicative factor associated with the local field enhancement of the nanocage structure (e.g., a hemisphere would introduce a factor of 3× via Schottky’s conjecture13,36,37) is not included for nanocages on the conical surface.

FIG. 1.

Typical arrays and emitters: (left) Conical array of emitters shown at edge-on perspective. Geometry of cones: height 1.6μm and base diameter 2.6μm; apex sphere designating 110 nm curvature (graphic not to scale) (after Fig. 1 of Ref. 53). (right) Top-down optical light image (courtesy of J. L. Shaw, NRL) of the emitter array for which side length of yellow triangle is 3μm and the diameter of the pink circle is 2.14μm.

FIG. 1.

Typical arrays and emitters: (left) Conical array of emitters shown at edge-on perspective. Geometry of cones: height 1.6μm and base diameter 2.6μm; apex sphere designating 110 nm curvature (graphic not to scale) (after Fig. 1 of Ref. 53). (right) Top-down optical light image (courtesy of J. L. Shaw, NRL) of the emitter array for which side length of yellow triangle is 3μm and the diameter of the pink circle is 2.14μm.

Close modal

Both triangular and square emitter grids or lattices will be treated below, with pitch given by multiples of the emitter height. The macroscopic shape of the array area will be taken as circular. For reference, a triangular grid with 3μm pitch (emitter density 1.28×107cm2) would require 78 nA per identical emitter to produce a 1A/cm2 current density. For typical temperatures of interest in heated cathodes (650°C900°C), the apex of the emitter may be field emitting, while the base regions may be thermally emitting. A general thermal-field model of current density will account for the transition between them.

Note that while the base of the true conical emitter tips (Fig. 1) makes an angle to the horizontal, their Point Charge Model (PCM) representation makes a right angle at the base (Fig. 2),resulting in vanishing fields there. Additionally, although the PCM is analytic, calculation of current over the conical emitter surface unavoidably requires numerical integration. To regain a rapid analytical model and the scaling relations it enables, the PCM method will be employed here to provide the parameters of an “effective” hyperbolic or ellipsoidal emitter matching the true conical emitter over much of the surface.

FIG. 2.

The point charges (left) and the resulting equipotential surface (right) shown in black for n=5 and r=0.7 for a pitch three times the height of an isolated emitter [that is, Δ=3Sn+1(r)=8.8235], and nt=2k+1=7. The radius of the base is approximately 1.61.

FIG. 2.

The point charges (left) and the resulting equipotential surface (right) shown in black for n=5 and r=0.7 for a pitch three times the height of an isolated emitter [that is, Δ=3Sn+1(r)=8.8235], and nt=2k+1=7. The radius of the base is approximately 1.61.

Close modal

In the PCM approach, a single emitter is composed of n point charges in a dipole arrangement37,54,55 where the charges are placed along the z^-axis at positions zj (j[1,n]), as in Fig. 2 (the companion negative charges below the z=0 plane are not shown). The background field E is embedded in Foq|E|=Uanode/D [compare Eq. (A4)], where q is the elementary charge, D is the anode-cathode separation, and the length scale lo is the emitter base radius corresponding to the n=1 dipole radius in the xy plane56 for values of a parameter r [defined in Eq. (1)] that governs the generation of the emitter, resulting in conical emitter shapes matching the zero equipotential of the electron potential energy U(x,y,z). In this convention, Fo is a “force due to the electric field” on an electron, but for brevity, it shall be referred to as “field.” It is convenient to render the calculations in a dimensionless format. If ρ~ and z~ are the dimensioned coordinates, then let ρ~=ρlo and z~=zlo, where ρ and z are understood to be dimensionless. Lastly, the locations of the point charges are placed such that

zj+1zjzjzj1r,
(1)

which implies that zj=zj1+rj and that the height of an isolated emitter is

j=0nrj=1rn+11rSn+1(r).
(2)

The factor r governs the sharpness and height of the emitter by dictating the relative magnitudes of the n point charges57 (see  Appendix B). The energy scale Folo is introduced. The dimensionless function Vn(ρ,z) is then defined by

Vn(ρ,z)=U(ρ~/lo,z~/lo)Folozj=1nλj(1rj+1rj),
(3)
rj±ρ2+(z±zj)2,
(4)

where U is the dimensioned potential energy, the radius factors rj± are distinct from the ratio factor r, and λj(r) is proportional to the magnitude of the charge at zj. λ absorbs factors such as q/4πε0 and others that otherwise would appear: it vanishes rapidly as j increases for small r, and, therefore, defining λj(r)rjPj(r) is numerically useful and defines Pj(r). A series of equations that uniquely determine Pj(r), particularly for the dipole case, as described previously,21,37,54,55 may be solved by solving the matrix relation of Eq. (B2).

Because Vn(ρ,z) is analytic, and the apex radius an and field enhancement βn exact [Eqs. (B4) and (B3), respectively], the total current from a single PCM emitter can be calculated via an integration over analytic functions. The field enhancement β is the ratio of the local field F on the surface [defined by zs(ρ)] to the background field Fo, and it is, therefore, dimensionless and may be calculated via gradients of Vn with z and ρ. Observe that an is dimensionless so that anlo is the dimensioned apex radius. It exponentially decreases with n for r<1 just as βn exponentially increases. Matrix methods for the evaluation of βn(r) and an(r) for the apex of an isolated dipole PCM emitter are summarized in  Appendix B. As r1, the shape of the resulting emitters approaches the shape of an ellipsoid: although r>1 solutions are possible, the line charge model58 becomes preferable in this limit. Conversely, for r<0.3, an isolated PCM emitter becomes disconnected from the base plane and begins to resemble a floating sphere albeit with a protruding bump at the apex; insofar as floating sphere models have relevance to field enhancement studies,59 a disconnected floating shape, even if unphysical, retains relevance as a means of modeling conditions at apex emission sites and relates to the utility of the prolate spheroidal models of  Appendix C.

The PCM with the parameters shown is motivated by but not dictated to conform with the physical structures of Fig. 1, as well as others described in the literature.53 Emphasis is placed on the analyticity of the model and its exactness with respect to the specification of the surface and the fields there. Greater correspondence to a physical emitter may be achieved through the addition of monopole charges at the base54 and abandonment of the self-similarity factor r relating the point charges, but that unnecessarily complicates the present model, which is designed to assess shielding, accuracy of work function characterization, and effects of microscale variation: the approximate correspondence is sufficient.

A “square lattice” is one for which emitters are placed on the intersections of a square grid so that the location of an emitter is characterized by its cartesian coordinates (xi,yj)=(x0+iΔ,y0+jΔ). Conversely, a “triangular lattice” is one for which emitters are at the vertices of equilateral triangles. The two arrangements are shown in Fig. 3 for emission sites within a (macroscale) circular boundary. In applications, the current density Jarray averaged over an array of emitters is sought, which is in turn related to an average current per emitter Itip and the array density of emitters ρarray such that Jarray=Itipρarray. For the same array density, a triangular (t) and square (s) array will have different pitches given by

Δs2=1/ρarray,Δt2=4/(3ρarray),
(5)

as shown by elementary arguments. Although the triangular array entails a larger pitch, more emitters are closer by, and so changes in pitch affect the degree of shielding as shown with line charge models.1,58 The location of the grid sites for a square lattice occupying a square total emission area is simply on a cartesian grid with grid elements of size Δ×Δ. In contrast, a triangular lattice is equivalent to two overlapping rectangular grids with grid elements of size Δ×3Δ with the origin of the first rectangular lattice at (0,0) and of the second at (Δ/2,3Δ/2).

FIG. 3.

Array of emitters designated by black open circles on a left square and right triangular lattice (S and T, respectively). Red circles in both cases bound a circular (C) area with a radius kΔ for k is an integer. Light gray lines show the lattice. Shown is a circular area for which the total emission area is within a circle (a square area would have all emitter sites with a square); the number of emitters along the x-coordinate axis is nt=2k+1, where k=3 for the arrays shown.

FIG. 3.

Array of emitters designated by black open circles on a left square and right triangular lattice (S and T, respectively). Red circles in both cases bound a circular (C) area with a radius kΔ for k is an integer. Light gray lines show the lattice. Shown is a circular area for which the total emission area is within a circle (a square area would have all emitter sites with a square); the number of emitters along the x-coordinate axis is nt=2k+1, where k=3 for the arrays shown.

Close modal

Consider a circular macroscopic emission area with the x-axis along the horizontal diameter. Let nt2k+1 be the number of emitters along that diameter, spaced every Δ (that is, xi=iΔ for kik). For a square (s) lattice, the placement of lattice sites above and below the diameter is governed by yj=jΔ for nt2i2jnt2i2 such that the diameter of the total emission area is ntΔ=(2k+1)Δ. For a triangular (t) lattice, a more general procedure (which can be generalized to treat an arbitrary number of randomly placed sites on a ring) is to introduce a matrix M^(ϕ) and a basis vector bij defined by

M^(ϕ)=[cos(ϕ)sin(ϕ)sin(ϕ)cos(ϕ)],bij=[δij0],
(6)

where the top and bottom components of bij are the x- and y-coordinates of a site, respectively, and the indices “do not” refer to matrix or vector elements. The factor δij is a radius defined by

δij={[jsin(π/3)]2+[ij(1cos(π/3))]2}1/2.
(7)

An initial rotation angle ϕij is defined by

tanϕij=jsin(π/3)ij[1cos(π/3)].
(8)

A hexagon of sites is then generated by initially rotating bij by M^(ϕij) followed by the application of M^(π/3) six times to find the sites of the six points of a hexagon. That is, the site of the lth vertex of the hexagon is located at M^(π/3)lM^(ϕij)bij. There are i hexagons for the ith ring (i.e., the first ring is composed of one hexagon giving 6 sites, the second of two hexagons giving 12 sites, and so on), where the rings are indicated by the red circles in Fig. 3. Such a procedure produces arrays with a hexagonal area, but for k6, these are the same as arrays allowing all sites within a circular region of diameter (2k+1). For k7, the procedure will begin neglecting emitters at the periphery, and so the designation “hexagonal“ (in contrast to “square”) shall be used in preference to “circular.”

Compared to an isolated emitter, the field enhancement β and apex radius a of an emitter in our array model changes depending on pitch Δ (an isolated emitter is equivalent to Δ). This is so because the hypothetical array is dictated by the zero-equipotential surface generated by the charge distribution, and the effects of distant charges on the field at a location contribute in a 1/r2 manner (where r here is the distance to the remote charge) and so are not necessarily negligible unless their distance is great. When the pitch decreases, the effects of these emitters on the shape of the equipotential increasingly contribute. Below, βn and an refer to the “apex” quantities for an isolated dipole PCM emitter composed of n point charges ( Appendix B); when the emitter is in an array, β(ρ) will refer to the local field enhancement along the surface of the rotationally symmetric emitter defined by zs(ρ), where ρ2=x2+y2, and a will denote the apex radius of the curvature for general pitch Δ. At the apex of the jth emitter, they are symbolically given by (generalizing Ref. 55)

β=zU(x,y,z)|z=zs,a=zU(x,y,z)ρ2U(x,y,z)|z=zs,
\vskip -8pt (9)

where U(x,y,z) is the (dimensionless) sum over all the contributing dipole PCM emitters on the lattice defining the array plus the background field and zs is the surface such that U(x,y,zs)=0. For an isolated dipole PCM emitter, U(x,y,z)Vn(ρ,z).

The presence of other PCM emitters complicates the determination of zs, which no longer conforms to the isolated emitter value z=Sn+1(r) due to the effect of adjacent emitters on the surface equipotential. Now, finding the surface requires using Newtonian iteration, with the initial estimate being Sn+1(r) for the center emitter on the axis. No loss of generality is entailed by focusing on those emitters along the x-axis; by symmetry, they are representative. The apex quantities will restrict consideration to the special case y=0, and focusing on the center emitter will ensure field enhancement is maximized at the apex and that the apex radius is unambiguous by symmetry, even though small departures from rotational symmetry occur.

A rapid scheme to find the surface makes use of second order finite difference approximations60 so that if fj=U(x,y,z+jδ) with δ1, then the kth iteration to find the surface is given by

zk+1=zk+2CB+B24AC,A=f12f0+f1,B=4δ(f1f1),C=8δ2f0,
\vskip 7pt (10)

with z0=Sn+1(r) at the apex or from a surface point found at a nearby location. The convergence is rapid, resulting in high accuracy by the eighth iteration. Thus, the determination of field enhancement along the x-coordinate is a two-step process: first, the surface is found by determining zs, and second, β(x,y,zs) is calculated. All the PCM emitters are created using identical charge distributions, but the equipotential surface is modified by which emitter on the lattice is under consideration: in what follows, the center emitter is always under consideration.

For the triangular lattice shown in Fig. 3, the profile of a typical center emitter in a PSS array was shown in Fig. 2, generated using n=5 and r=7. The surface equipotential generated by the above iterative evaluation is shown in Fig. 4 for varying numbers of emitters along the x-axis nt labeled by k, where nt=2k+1 as in Fig. 3. Off-center emitters or emitters on a square lattice are visually similar but not identical to the center emitter due to altered shielding, in a manner suggested by Fig. 4.

FIG. 4.

Surface equipotential for increasing array size. Number of emitters along the x-axis for a triangular lattice for a hexagonal array is nt=2k+1. All other parameters are the same as those used for Fig. 2.

FIG. 4.

Surface equipotential for increasing array size. Number of emitters along the x-axis for a triangular lattice for a hexagonal array is nt=2k+1. All other parameters are the same as those used for Fig. 2.

Close modal

The PCM is used to parameterize a prolate spheroidal model to enable the usage of analytic I(Fo) relations and suppress negligible variations at the surface induced by a discrete array of distant emitters. Compared to a floating sphere as in Fig. 5, an ellipsoidal and hyperbolic shape compare well to the PCM surface with a radius of curvature specified by Eq. (9): the spherical emitter is a poor approximation to the body of the emitter. With respect to floating sphere or simple analytic models relying only on apex height and tip radius, the remainder of the emitter structure (and the base on which it rests) unavoidably plays a role in defining the overall field enhancement,61 as shall be seen. The prolate spheroidal models11,33–35 perform better and enable analytic field enhancement factors,23,24 but the boundary conditions (e.g., anode-cathode gap or isolated emitter) do not generalize to arbitrarily spaced anodes or an array of similar emitters.

FIG. 5.

The surface equipotential (dots) for n=5, r=0.7, and k=3 compared to the spherical (green), ellipsoidal (blue), and hyperbolic (red) fits based on Eqs. (11) and (12).

FIG. 5.

The surface equipotential (dots) for n=5, r=0.7, and k=3 compared to the spherical (green), ellipsoidal (blue), and hyperbolic (red) fits based on Eqs. (11) and (12).

Close modal

The prolate spheroidal representation becomes useful when the PCM is used to parameterize it. The parameters associated with an ellipsoidal (ell) and hyperbolic (hyp) emitter, as well as the spherical (sph) one, are obtained. Let the numerically found equipotential surface coordinates be characterized by the pairs [xi,zs(xi)] (y=0 for convenience so that ρx), let z1=zs(0) and z2=zs(x2). The pairs [xj,zj=zs(xj)] define the numerically determined points of the surface and x10 corresponds to the apex coordinate. In terms of them, and based on the relation z(x) for each of the coordinate systems as in  Appendix C, define

asph=z1z22+x222(z1z2),aell=x2z1z12z22,ahyp=(z1z2)(3z1z2),
\vskip 4pt(11)

in terms of which

zs(x)|sph=z1[asphasph2x2],zs(x)|ell=z1aellaell2x2,zs(x)|hyp=2z1ahypx2x2+(x2z1/ahyp)2.
(12)

Of the three representations, the hyperbolic representation provides the best correspondence with the physical C12A7-on-PSS emitters by softening the sharp (90°) angle that the emitter makes with the base plane. The surface fields are better represented by the hyperbolic parametrization.

The field enhancement factor β(x) as a function of position x along the surface for a hyperbolic emitter is analytic (ρx for y0). That it is a reasonable approximation to the PCM β(x) is demonstrated as follows. The hyperbolic emitter field enhancement factor62 is given by

β(x)β(0)|hyp=sinυosinh2η+sin2υo,
(13)

where υo specifies the shape of the emitter and the coordinates are as in  Appendix C. Coupled with the hyperbolic representation derived from Eq. (C2) that x=sinhηsinυo, then

sinυo=(x2β2/β1)1/2[1(β2/β1)2]1/4,
(14)

where βj=β(xj) and x1=0 as before. Using Eq. (14) in Eq. (13) shows that for the (n,r) parameters chosen, a good overlap between the hyperbolic model and the PCM array model is obtained as a consequence of the closeness of the zero equipotential of Fig. 2 to a hyperbolic shape for the range 0.65r0.85 for n5. If the hyperbolic model is to be a good approximation, then n and r must be chosen appropriately together: [(n,r)=(5,0.7) is such a choice, as is (n,r)=(12,0.85)]; other combinations such as (n,r)=(3,0.95) are better represented by the ellipsoidal model which would necessitate deriving the appropriate replacements to Eqs. (13) and (14), e.g.,

β(x)β(0)|ell=sinhηocosυsinh2ηo+sin2υ,
(15)

with sinhηo determined analogously to Eq. (14). The hyperbolic approximation is the better choice for comparisons to the C12A7 emitters.

The finite area of emission associated with a unit cell (the yellow hexagonal region containing one emitter, as in Fig. 6) and the shielding are affected by the pitch factor Δ and the size of the macroscopic array. Although the PCM method is effective for small values of nt (e.g., k ranging from 0 to 10), a circular emission area with a diameter of 1 cm and a pitch of Δ=5μm entails k1000. Emitters around the rim have negligible area of contribution: the outer edge of emitters for a D=1 cm diameter disk accounts for (Δ/D)(2Δ/D)<0.01% of the emission area. A large area array can, therefore, be taken as the k limit with respect to shielding effects. The βk limit may be estimated by observing that βk/β0 is to a good approximation linear in 1/k, as the number of emitters in each successive ring scales with k, while the contribution to the field declines faster than 1/k2 for dipoles. The calculations behind Fig. 7 for Δ=3Sn+1(r) for r=0.7 and n=5 result in the p=3 line of Table I. Performing the same analysis with increasing pitches, indexed by p=Δ/Sn+1(r) for n=5 and r=0.7 [for which S6(0.7)=2.9412], shows the diminishing effect of shielding as Δ increases. As a result, PCM models capture how asymptotic shielding effects scale with k and, therefore, justifies use of the PCM as a design tool for pitch spacing. For example, it is seen that the enhancement for well-spaced arrays (p6) is within 4% of the unshielded case (close to the p=12 value).

FIG. 6.

Areas associated with emission regimes for a pitch of Δ and emitter base radius of δ. For sufficiently high fields, the blue center circle is the area of field emission, the green center circle the area of thermal-field emission, the orange circle the area of thermal emission from the tip, and the yellow hexagon the area of thermal emission from the base. From the side, the emitters would appear to have a hyperbolic shape.

FIG. 6.

Areas associated with emission regimes for a pitch of Δ and emitter base radius of δ. For sufficiently high fields, the blue center circle is the area of field emission, the green center circle the area of thermal-field emission, the orange circle the area of thermal emission from the tip, and the yellow hexagon the area of thermal emission from the base. From the side, the emitters would appear to have a hyperbolic shape.

Close modal
FIG. 7.

Field enhancement factor βk in arrays characterized by nt=2k+1 compared to an isolated emitter β0 (k=0) for a pitch set equal to Δ=3Sn+1(r).

FIG. 7.

Field enhancement factor βk in arrays characterized by nt=2k+1 compared to an isolated emitter β0 (k=0) for a pitch set equal to Δ=3Sn+1(r).

Close modal
TABLE I.

The asymptotic values of βk/β0 as k for n=5 and r=0.7 for values of pitch p=Δ/Sn+1(r), with Sq=squarearea, Hex=hexagonalarea, Car=cartesianlattice, and Tri=triangularlattice.

pSq CarHex CarSq TriHex Tri
0.309 553 0.305 726 0.259 846 0.259 088 
0.559 134 0.557 655 0.526 304 0.525 689 
0.657 410 0.656 726 0.639 896 0.639 570 
0.698 481 0.698 116 0.688 609 0.688 429 
0.717 969 0.717 752 0.711 976 0.711 869 
0.728 255 0.728 117 0.724 379 0.724 310 
0.734 156 0.734 063 0.731 516 0.731 470 
0.737 770 0.737 705 0.735 897 0.735 864 
10 0.740 103 0.740 055 0.738 727 0.738 703 
11 0.741 674 0.741 637 0.740 634 0.740 616 
12 0.742 769 0.742 741 0.741 965 0.741 951 
pSq CarHex CarSq TriHex Tri
0.309 553 0.305 726 0.259 846 0.259 088 
0.559 134 0.557 655 0.526 304 0.525 689 
0.657 410 0.656 726 0.639 896 0.639 570 
0.698 481 0.698 116 0.688 609 0.688 429 
0.717 969 0.717 752 0.711 976 0.711 869 
0.728 255 0.728 117 0.724 379 0.724 310 
0.734 156 0.734 063 0.731 516 0.731 470 
0.737 770 0.737 705 0.735 897 0.735 864 
10 0.740 103 0.740 055 0.738 727 0.738 703 
11 0.741 674 0.741 637 0.740 634 0.740 616 
12 0.742 769 0.742 741 0.741 965 0.741 951 

The present analysis did not accommodate changes in the apex radius an(r) with pitch Δ (in contrast to prior studies using the LCM1,29,58 where the apex radius was held fixed): as the equipotential curves become progressively less blunt as Δ decreases, βk (which is proportional to the inverse apex radius) correspondingly decreases. As an example, for Δ=1.85Sn+1(r) and k=3, the base of the dipole PCM emitters is touching in a way that is reasonably approximated by cos2(x/2) scaled to the height of the emitters (and, therefore, approximately represents roughness models using trigonometric functions, e.g., Ref. 20). Without adjusting conditions so that weak variation in the apex radius can result, the relationship between the apex radius an(r) and pitch is shown in Fig. 8 for n=5 and k=3 and demonstrates that for sufficiently large pitch, the changes to the apex radius are weak. If instead apex radius were to be engineered to remain constant, the decline in βk with decreasing Δ would be less.

FIG. 8.

Apex radius factor asph(r) in arrays characterized by n=5 and k=3 as a function of p=Δ/Sn+1(r) (ratio of pitch to isolated emitter height).

FIG. 8.

Apex radius factor asph(r) in arrays characterized by n=5 and k=3 as a function of p=Δ/Sn+1(r) (ratio of pitch to isolated emitter height).

Close modal

Because several factors determine the apex radius, care in how they are altered is required. For example, increasing the apex point charge may appear preferable to altering r, as changes to the apex radius can be made by changing the magnitude of the topmost point charge without affecting the emitter shape at the base. Studies that preserve the apex radius (but endure an enlarging of the base radius as a consequence) using a line charge model (LCM)1,24,29,58 anticipate what may be expected for a small array of PCM emitters. Changes to the topmost point charge to preserve radius will be the subject of a separate study, but here, the relation of the PCM charges λj is held fixed to those of an isolated emitter to allow for unambiguous quantitative comparisons.

Such a distinction reflects the importance of macroscale, mesoscale, and microscale surface features.1 Altering the topmost point charge effectively mimics changing the microscale protrusions without (greatly) altering the mesoscale base radius, whereas altering r (or altering the line charge length in LCM representations of elliptical emitters) effectively mimics changing the mesoscale dimensions without altering the microscale protrusions if apex curvature is held fixed. In the analogous LCM models, changes in the line charges simultaneously affect macroscale and mesoscale features such that holding a parameter like the apex radius fixed requires monitoring.

The behavior of β(x) in Fig. 9 anticipates severe problems in modeling total current from an array using only the canonical emission equations,24 namely, the FN equation (field emission) and the RLD equation (thermal emission), when work function and field enhancement are such that the apex may experience field emission when the flat regions between emitters are dominated by thermal emission. These problems are grouped as follows:

  1. when both field and thermal emission contributions are present, then for a portion of the surface (the side walls of the hyperbola), thermal-field emission contributes and is not well-modeled by FN and RLD or their sum;

  2. when the pitch Δ is small, then the field at the surface of the flat regions between emitters is “smaller than the background field,” and, therefore, inclusion of the Schottky factor (ΦΦ4QFo in the notation of Ref. 24) in the RLD equation overestimates the current;

  3. inferring Φ from thermal data using RLD plots or the field data using FN plots is mistaken because, first, the Schottky factor is overestimated between emitters by neglecting shielding effects in the former, and the apex field enhancement factor is reduced by shielding due to adjacent emitters in the later; and lastly,

  4. the thermal-field (TF) current density substantially exceeds both JRLD and JFN in the transition region where thermal and field effects are comparable.

FIG. 9.

Field enhancement factor β(x) along the surface shown in Fig. 5. The hyperbolic line βh(x) shows the field enhancement factor evaluated using the approximation of Eq. (13) along the hyperbolic surface.

FIG. 9.

Field enhancement factor β(x) along the surface shown in Fig. 5. The hyperbolic line βh(x) shows the field enhancement factor evaluated using the approximation of Eq. (13) along the hyperbolic surface.

Close modal

A hyperbolic representation parameterized by the PCM method becomes the preferred method for analyzing thermal-field and shielding effects because it (i) enables a highly accurate analytical approach, (ii) allows a clean demarcation of the regions where thermal, thermal-field, and field emission dominate, (iii) provides a semianalytic means of ascertaining when current-voltage data transition from one regime to another,63 and (iv) provides an accurate means of assessing the consequences of various configurations on the estimates of work function and emitter parameters that result from characterization efforts. The numerical integration of a general thermal-field current density JGTF over the surface of a PCM emitter will be compared to the thermal, thermal-field, and field contributions of the PCM-hyperbolic model in which an emitter of base radius δ is situated at the center of a circle of diameter Δ. The current density is

JGTF(F,T)ARLDT2N(n,s),n(F,T)βTβF(Em),s(F,T)βF(Em)(Eoμ),
(16)

where the terms and their evaluation are described in  Appendix A. The equation becomes, in the appropriate limits, the RLD equation of thermal emission in Eq. (A20) and the FN equation of field emission in Eq. (A21). Four regions are characterized by limits of JGTF and areas suggested in Fig. 6 if conditions are such that the onset of field emission occurs at the apex and are denoted field (F), thermal-field (TF), thermal (T) on the emitter, and thermal (To) on the base, as labeled in Fig. 6:

  • Current density at the apex is characterized by n>1 for which if F is high, JGTFJFN.

  • A region adjacent to the apex is characterized by n=1 for which N(1,s)=(s+1)es.

  • Current density near the base of the emitters is characterized by n<1 for which JGTFJRLD.

  • Hexagonal base (flat) region outside the thermal emitter region for which the emission area is the difference between the area of a hexagon and a circle or (33/2)Δ2πδ2.

For each of these conditions, the numerical evaluation of JGTF(F,T) is performed by Eq. (A16) when n=βT/βF1 and Eq. (A17) when n=1. Close inspection, though, reveals that the n1 limit of Eq. (A19) differs from Eq. (A17), as can be seen by writing the former as

N(n,s)=[No(n,s)+ΔN(n,s)]es,No(n,s)=(e(n1)s1)Σ(n),ΔN(n,s)=n2Σ(1n)+Σ(n).
\vskip -4pt(17)

At n=1, No(1,s)=s exactly because the vanishing of e(n1)s1 offsets the singularity in Σ(n) as n1. It is seen that

limn1ΔN(n,s)=7π4180+π236=1.0780,
(18)

which exceeds unity. The accuracy of ΔN, therefore, determines the accuracy of Eq. (A19). However, s is large: when n=1, then s40 for metal-like work functions (Φ=4.5eV) and s10 for low work functions (Φ=2eV) when F=1eV/nm. As a result, the small error (ΔN1)/s1 entailed by Eq. (18) is not easily discerned visually. Therefore, the special usage of Eq. (A17) directly when n=1 instead of Eq. (A19) as n1 will not result in excessive visual incongruities when using the prescriptions of  Appendix A with respect to the relation of T to Tmax and Tmin in the evaluation of n and s.

Consideration of the rotationally symmetric central emitter next requires the evaluation of differential surface elements in the shape of circular ribbons of radius ρj for the jth ribbon. There is no requirement that ρj be evenly spaced: near the apex of the emitter where fields vary most strongly, ρj is more densely spaced than further from the axis of symmetry. Field enhancement values βj=β(ρj) may be constructed from Eq. (15) with the replacement xρ. The current density at the locations ρj is then defined by the components

Jj=JGTF[β(ρj)Fo,T].
(19)

The area of a circular ribbon of diameter ρ and thickness dl is dA=2πρdl. For a surface defined by zs(ρ), where dl2=dρ2+dz2, then

dA=2πρdρ2+dz2=2πρdρ1+(ρzs)2.
(20)

For a surface zs(ρ) specified by the discrete points (ρj,zj) with zj=zs(ρj), and approximating the gradient of zs(ρ) using finite differences,60 then

ΔAj+1/2πlo2(ρj+12ρj2)[1+(zj+1zjρj+1ρj)2]1/2,
(21)

where zj=zs(ρj). The emitted current from the differential surface element is then

dIj+1/2=12[JGTF(ρj)+JGTF(ρj+1)]ΔAj+1/2,
(22)

where the notation JGTF(ρ) is the current density at (ρ,zs(ρ)) experiencing a β-factor of β(ρ). The total current within a circle of radius ρN+1 is then

Ipcm=j=1NdIj+1/2.
(23)

The usage of the hyperbolic surface zh(ρ) simplifies the model but introduces an ambiguity indicated in Fig. 9: at some point, βh<β(x), whereas the intermediate region between cones experiences a surface field less than but comparable to Fo, and current density well-characterized by thermal emission will be affected by Schottky barrier lowering for the flat regions past the edge of the hyperbolic cone shown in Fig. 5 corresponding to the yellow hexagon of Fig. 6. For simplicity, the current density in the yellow region will be evaluated for β(ρ)β(Δ/2) (the field enhancement of the midpoint between cones) as an upper bound approximation. The current from the unit cell (the hexagon of Fig. 6) is then the sum of Ipcm and the thermal current from the hexagonal region or

IunitIpcm+J[β(Δ/2)Fo,T](332Δ2πδ2),
(24)

where for purposes of estimating the emission parameters, the equality will be assumed.

Simulations were performed using parameters for which the PCM unit emitter was well approximated by a hyperbolic geometry. Such parameters give rise to emitters and pitches similar to those in Table II. The simulation parameters were designed to bring out differences between GTF, FN, and RLD current density estimates for a mesoscopic area unit cell with microscopic field emission sites using a hyperbolic representation of the emitter. Higher accuracy is achievable for the emitters of Fig. 1, but that would require more intensive computations with respect to the PCM and usage of actual β(x) and zs(x) as determined numerically and so shall be part of a separate study. The present approach is sufficient to assess the accuracy of predicting the emitter parameters in light of transitions through the thermal-field regime along unit cell emitters. A moderate tip number parameter ntip=2k+1 of k=3 was chosen for computational speed while still providing representative results.

TABLE II.

Values and parameters used in the characterization of current-voltage relations for conical emitters subject to thermal-field emission regions.

ParameterSymbolValueUnit
Temperature T 1173 
Work function Φ 2.3 eV 
Fermi energy μ eV 
Length unit lo μm 
Charge number n  
Self-sim ratio r 0.7  
Pitch Δ 8.8235  
Tip number k  
Height zs(0) 2.9412  
Base radius δ 1.7071  
Lattice  Tri  
Area  Hex  
ParameterSymbolValueUnit
Temperature T 1173 
Work function Φ 2.3 eV 
Fermi energy μ eV 
Length unit lo μm 
Charge number n  
Self-sim ratio r 0.7  
Pitch Δ 8.8235  
Tip number k  
Height zs(0) 2.9412  
Base radius δ 1.7071  
Lattice  Tri  
Area  Hex  

The total emitted current Iunit from a unit cell (one emitter at the center of a hexagonal area as in Fig. 6) was evaluated using Eq. (24) with J=JGTF(F,T). The discrete factors xj=(j/Nx)2(Δ/2) with Nx=128 so as to more finely grid near the apex of the emitter. The background field Fo is taken to be linearly proportional to an anode potential φa (e.g., Fo=qφa/D for which φa=2kV if D=100μm and Fo=0.02eV/nm). Comparisons were made to standard models using the Fowler-Nordheim equation [Eq. (A21)] for emission from the apex over a notional emission area62 of πas2g[β(0)Fo] where for hyperbolic emitters,23,64

g(F)=cos2θb+νsin2θ,
(25)

where65 b(F,Φ)=42mΦ3/(3F) and ν(Φ)=(8Q/9)2m/Φ. For a cone angle of tan(θ)=1.3/1.6 (or θ=39°) and Φ=2.3eV, g(F) simplifies to g(F)=0.88150F/(34.869F) for F in units of eV/nm. Using the RLD equation [Eq. (A20)] instead of JGTF, then for the same total area characterized by Φ and a surface field of β(Δ/2)Fo where β(Δ/2)=0.7121 for the parameters of Table II, the IRLD result is obtained. Using the FN equation [Eq. (A21)] instead of JGTF for a notional emission area governed by Eq. (25) with F=β(0)Fo, then the IFN result is obtained. Comparisons are made using a simple I vs F plot, the same plot in RLD coordinates [ln(I) vs F], and in FN coordinates [ln(I/F2) vs 1/F], as shown in Figs. 10, 11, and 12, respectively. Departures of the GTF (points) from linearity in either FN coordinates or RLD coordinates (lines) indicate the onset of thermal-field effects.

FIG. 10.

Current Iunit (nA) from a unit cell as a function of background field Fo (eV/nm) (symbols) for parameters of Table II. Lines correspond to assuming all current is from field emission from the apex over a notional area 2πasph2g(βFo)JFN(βFo), with β=β(0) being the field enhancement at the apex.

FIG. 10.

Current Iunit (nA) from a unit cell as a function of background field Fo (eV/nm) (symbols) for parameters of Table II. Lines correspond to assuming all current is from field emission from the apex over a notional area 2πasph2g(βFo)JFN(βFo), with β=β(0) being the field enhancement at the apex.

Close modal
FIG. 11.

Same as Fig. 10 but in FN-related coordinates, wherein ln(Iunit/Fo2) is linear in 1/Fo with Fo proportional to the anode potential.

FIG. 11.

Same as Fig. 10 but in FN-related coordinates, wherein ln(Iunit/Fo2) is linear in 1/Fo with Fo proportional to the anode potential.

Close modal
FIG. 12.

Same as Fig. 10 but in RLD-related coordinates, wherein ln(Iunit) is linear in Fo.

FIG. 12.

Same as Fig. 10 but in RLD-related coordinates, wherein ln(Iunit) is linear in Fo.

Close modal

Consider now the reverse problem of determining emitter parameters from the Iunit(V) values in a manner comparable to the usage of experimental data to determine work function. The effective work function Φeff is determined from the slope Φeff/kB on an RLD plot of ln(Iunit/T2) vs 1/T. Such an analysis is shown for Table II parameters in Fig. 13 (symbols). Were JGTF(F,T) replaced by JRLD(Fo,T) with Φ=2.3eV, the gray lines would result and the slopes would correspond to ϕ=Φ4QFo. With the field enhancement at the apex of β=7.6202, the apex of the emitter increasingly contributes field emission current, and the I(Fo) data demonstrate how that contribution reduces the slope and, therefore, underestimates the fitted work function.

FIG. 13.

Estimation of Φeff from the slope of Iunit on an RLD plot for Fo=(0.02,0.06,0.10)eV/nm gives Φeff=(2.2959,2.1398,1.7551)eV, respectively. Units of Iunit is nA and T in K. All “rld” slopes (gray lines) correspond to Φ=2.3eV.

FIG. 13.

Estimation of Φeff from the slope of Iunit on an RLD plot for Fo=(0.02,0.06,0.10)eV/nm gives Φeff=(2.2959,2.1398,1.7551)eV, respectively. Units of Iunit is nA and T in K. All “rld” slopes (gray lines) correspond to Φ=2.3eV.

Close modal

For field data, the effective field enhancement factor βeff is determined (if Φ is known) from the slope on an FN plot of ln(Iunit/Fo2) vs 1/Fo (even though there are problems of interpretation in doing so,66,67 and factors such as resistance along the emitter are ignored but are experimentally relevant). The slope on an FN plot is now approximately given by b(Fo,Φ)/β, where b(F,Φ) is as encountered in Eq. (25) and the variation of the other factors in IFN with β is weak and, therefore, ignored. Such an analysis is shown for Table II parameters in Fig. 14 (symbols). Were JGTF(F,T) replaced by JFN(βFo), the single gray line would result and the slope would correspond to β=7.6202. The sides of the emitter and the hexagonal region of the base increasingly contribute thermal emission current, and the I(Fo,T) data demonstrate how that contribution reduces the slope and, therefore, underestimates the fitted field enhancement factor.

FIG. 14.

Estimation of βeff from the slope of Iunit on an FN plot for T=(77,300,600)K gives βeff=(7.6670,7.8991,10.667), respectively. Units of Iunit is nA and Fo is eV/nm. Single “fn” line (gray) corresponds to β=7.6202.

FIG. 14.

Estimation of βeff from the slope of Iunit on an FN plot for T=(77,300,600)K gives βeff=(7.6670,7.8991,10.667), respectively. Units of Iunit is nA and Fo is eV/nm. Single “fn” line (gray) corresponds to β=7.6202.

Close modal

The reverse analysis in which slopes on an RLD and FN plot are used to infer work function and field enhancement, respectively, indicates quite clearly that if emission conditions are such that a thermal-field contribution to the emitted current is present, then the procedure is undermined. Furthermore, when the emitter contains surface structure that introduces a microscopic contribution to a mesoscale emission area, and when the space between the C12A7 cones under the present study composed of squat and close-packed cylindrical emitter shapes on which nanocage structures occur,52 then the methods by which either the work function or the field enhancement factor is estimated and total emission current deduced theoretically matters significantly in both the characterization of such cathodes and predictions of their performance.

The random surface roughness model will be motivated by photocathode surfaces. Existing photocathodes for next generation and x-ray Free Electron Lasers (xFELs) have difficulty meeting future requirements, motivating efforts to understand cathode intrinsic emittance and limitations to performance.10,16,19,20,41,43,68,69 Predictions of emittance require estimating the MTE and, therefore, a model of trajectories close to the cathode surface. Two models are developed for that purpose. In the first, a Surface Roughness Model (SRM) based on a modified point charge model is developed to evaluate required analytical β-factors and differential surface elements. In the second, a ballistic-impulse model is developed for the MTE estimates. The need for both is because intrinsic emittance (that due to the photocathode itself) is emerging as an irreducible and a significant component of the emittance of an electron beam, and the random nature of surface roughness creates difficulties for analytical models.

The PCM introduced in Sec. II requires some modifications in order to apply to a surface roughness model, chiefly being that the locations of the point charges are now randomly placed with random magnitudes (compare to single bump and sinusoidal treatments in Refs. 10, 20, and 69), and the undulations are less pronounced than for the conical emitters considered in fabricated arrays of Sec. II. Now, a placement of Nx2 point charges on a grid is considered. On a regular (nonrandom) array, the (xi,yj) coordinates of the (i,j)th point charge, using the convention that (x,y,z) are scaled by L and so are dimensionless, are at

xi(2iNx1)2(Nx1)(12δ),yj(2jNx1)2(Nx1)(12δ),
(26)

so that the simulation region is centered at (x,y)=(0,0) and spans 1/2(x,y)1/2, and the outermost point charges are a distance δ from the boundary, as in Fig. 15. In contrast to how the PCM is used to model conical emitters as in Sec. II, here, only a single point charge of varying magnitude qij is used, placed at z=zo, with its image charge at zo (dipole model). To randomize the locations, a radius and an angle characterized by η1 and 2πη2 are added such that xixi+η1cos(2πη2) and yjyj+η1sin(2πη2), where the ηj are uniformly distributed random numbers 0ηj1. The potential energy is then

U(r)=FozL+q4πϵ0Li,jqij{1rij1rij+},
(27)
(rij±)2=[(xxi)2+(yyj)2+(z±zo)2].
(28)
FIG. 15.

Location of qij without (green ) and with (red ) random offset for Nx=4: see Eq. (26). δ=0.125.

FIG. 15.

Location of qij without (green ) and with (red ) random offset for Nx=4: see Eq. (26). δ=0.125.

Close modal

To closer approximate a large surface, periodic boundary conditions are used: that is, the pattern of point charges shall be periodically repeated every L, or xNx+1=x1 and similarly with yj. Analogous to Eq. (3), the dimensionless equation for V(x,y,z) becomes

U(r)FoL=z+i,j,k,lλij{1rijkl1rijkl+},
(29)

where

(rijkl±)2=(xxi+k)2+(yyj+l)2+(z±zo)2,
(30)

where U is the (dimensioned) potential energy, and the summations over k and l range from n to n [as there is only a single charge and its image now, the variable n is now a measure of the number of repeating tiles so that the area of the surface is (2n+1)2L2]. In practice, n=2 is sufficient to mimic infinite lattice periodicity when zo is small.

The surface is found by locating the points for which V(x,y,zs)=0, which defines zs(x,y). Newtonian iteration is used to find zs via a finite difference representation of zV[V(z+ϵ)V(zϵ)]/2ϵ so that the equation V(z+ϵ)=0=V(z)+ϵzV is iteratively solved by

zszs2ϵV(zs)V(zs+ϵ)V(zsϵ),
(31)

where the dependence of V on x and y is suppressed. The initial (first guess) value of zs is where V is approximately linear in z. Although the single dipole point charge model allows for the gradients to be calculated analytically, the finite difference method is sufficiently accurate and is moreover rapid and, consequently, useful when the evaluation of O(104) values of zs is sought, as in Fig. 16.

FIG. 16.

The surface zs(x,y) as found using Eq. (31) (zs is dimensionless so the color scale is unitless).

FIG. 16.

The surface zs(x,y) as found using Eq. (31) (zs is dimensionless so the color scale is unitless).

Close modal

The field F(x,y,z) is also analytic, but as before, the evaluation of the components of F=(Fx,Fy,Fz) using finite difference methods is rapid. Therefore,

2ϵFjFoV(xj+ϵ)V(xjϵ),
(32)

where xj is x, y, or zs(x,y). The total force F2=Fx2+Fy2+Fz2 is the product of a field enhancement factor β(x,y,zs) with the background field Fo such that |F|=βFo. The variation of β(x,y,zs) for the surface of Fig. 16 is shown in Fig. 17, and it is used to find the local field |F|=βFo used in the emission equation JGTFP(F,T). The enhancement is seen to be above background (β>1) near the apexes of the protrusions, but significantly, “suppression” (β<1) occurs in the valleys between the protrusions and results in reduced emission compared to a planar surface. The nature of field suppression will be taken up separately.70 

FIG. 17.

Field enhancement factor β(x,y,zs) using Eq. (32) (β is dimensionless so the color scale is unitless).

FIG. 17.

Field enhancement factor β(x,y,zs) using Eq. (32) (β is dimensionless so the color scale is unitless).

Close modal

Curvature of the random surface entails that the differential surface element dS is larger than the product δxδy. Define the area factor A(x,y)=δS/(δxδy), which is related to the gradients of zs(x,y) by (assuming δx=δy)

A(x,y)=dS(δx)2=1+(xzs)2+(yzs)2,
(33)

where finite difference approximations may be used to evaluate the gradient terms (Fig. 18). Because in simulation, the small but finite sized grid squares approximate the differential surface elements, in calculation, Aij takes the place of A(x,y). The product of the current density J on that grid element with A then results in the current element I(x,y)=J(x,y)A(x,y)δx2. When performing a trajectory analysis, it will be more convenient to weight the trajectories by a number w(x,y) which gives the number of electrons emitted from A(x,y) compared to the number emitted from the minimum current element or

w(x,y)wmin=I(x,y)IminJ(x,y)A(x,y)Jminδx2.
(34)
FIG. 18.

The area element factor A(x,y) using Eq. (33). [A(x,y) is dimensionless so the color scale is unitless].

FIG. 18.

The area element factor A(x,y) using Eq. (33). [A(x,y) is dimensionless so the color scale is unitless].

Close modal

For photoemission, if J(x,y) is evaluated using the Fowler-DuBridge limit of Eq. (A16) (where n21),

J12(ωϕ)2+(πkBT)26,
(35)

on each element with Φ=2.3 eV, T=300 K and Fo=0.05 eV/nm, 4QFo=0.268325 eV, and λ=532 nm (ω=2.33053 eV), then Fig. 19 results. The “hole” in the center is due to the random placement of the four center emitters of Fig. 15 all being away from the center and is, therefore, a random effect. Because field suppression is most pronounced in that hole, the difference between the local work function and photon energy is narrowed further, resulting in a large ratio of w/wmin near the apexes.

FIG. 19.

Weighting factor w(x,y) for JJFD with Φ=2.3eV, T=300K, and λ=532nm. The number of electrons emitted from each area element A(x,y) is proportional to w(x,y).

FIG. 19.

Weighting factor w(x,y) for JJFD with Φ=2.3eV, T=300K, and λ=532nm. The number of electrons emitted from each area element A(x,y) is proportional to w(x,y).

Close modal

The Ballistic-Impulse model ( Appendix D) enables a specification of velocity components at the anode plane acting as the virtual cathode in a unit cell model for beam optics codes such as MICHELLE. If that plane is at z=za, then p=pa is the time taken for the trajectory to get to the plane and is evaluated from Eq. (D13) as the solution of a quadratic equation in p, from which vρ and vz at ta=pa(a/vo) can be evaluated.

A sufficiently planar equipotential above the surface in Fig. 16, functioning as an “anode plane,” may be treated as a virtual particle source for beam optics codes if the distribution of electrons in velocity and position can be found. That in turn requires the trajectories of electrons as they cross that plane to be well approximated, and the ballistic-impulse approximation is a means to do so. It first finds an adequately planar equipotential of the potential energy U and then launches electrons from the surface zs(x,y) with an impulse-modified initial velocity, after which the electrons follow ballistic trajectories in a constant field Fo.

A measure of the flatness of an equipotential surface is given by the variation of departures from the mean value. Let the highest point of the surface be zc=max(zs(x,y)), and let za be the intended location of the virtual anode plane. Define zj=za(zazc)(j/N) as intermediate locations: the variation of the rms departure as a function of j is then given by

ΔVj=Vj2Vj21/2,
(36)

where VjV(x,y,zj) and the mean Vj is approximated by a sum over uniformly distributed values over the surface (a grid). As shown in Fig. 20, undulations reduce to approximately 1% of Vj by z2zc. Consequently, the anode plane that acts as the virtual cathode can be taken at a height above the surface corresponding to the length scale dictated by the surface roughness itself.

FIG. 20.

The decline of ΔVj is 1.28% of Vj for z2zc, where zc is the highest point on the surface. The fall-off is well approximated by yxp with p>4.

FIG. 20.

The decline of ΔVj is 1.28% of Vj for z2zc, where zc is the highest point on the surface. The fall-off is well approximated by yxp with p>4.

Close modal

Next, the evaluation of β(x,y) on the surface zs(x,y) is found in order to construct the impulse factor σ used to alter the initial velocity of vo, assumed normal to the surface and, therefore, parallel with the field line. The local impulse factor, however, requires specification. The observation that β(θ)=3cosθ for the hemisphere suggests that the generalization to Eq. (D15) should be

σ(x,y)1+12Cκβ(x,y).
(37)

The analytic approximation of (vx,vy,vz) (in units of |vo|) at z=za is then specified by the sequence

  1. Find zs(x,y) using the PCM via Eq. (31).

  2. Evaluate β(x,y) at zs.

  3. Specify σ and κ via Eqs. (37) and (D8), and C by comparisons to trajectories [or Eq. (D18)].

  4. Evaluate v(t) for launches distributed over surface.

  5. Specify the time ta, defined by z(ta)=za2zc.

  6. Weight v(ta) by w(x,y)/wmin via Eq. (34) and construct v-distributions.

Regarding Item 4, the components vx and vy are unchanged from their launch values in the Ballistic-Impulse approximation: vk(x,y)σ(x,y)vok^ for the k^th component, which allows vz(ta) in the Ballistic-Impulse approximation to be given by Eq. (D17), rewritten as

vz(ta)2(σvz(0))2=2Fom[zazs].
(38)

Because the mesoscale L is already small, smaller errors in (x,y) initial positions are ignorable for the macro scale beam dynamics. Tracking only the velocity components without retaining x(ta) or y(ta) is, therefore, permissible and simplifies computational demands. If the surface normal is n^ such that n^z^cosθ=Fz/βFo, and so, vz(0)=voFz/βFo=vocosθ where71FzL=zU|z=zs, then vx(0)=vosinθcosϕ and similarly for vy, to give

vz(0)vo={2κ(zazs)za+(σcosθ)2}1/2,vx(0)vo=σsinθcosϕ,vy(0)vo=σsinθsinϕ,
\vskip 4.6pt (39)

where ϕ, like θ, is defined with respect to the surface normal n^ and all initial quantities are defined at (x,y,zs).

Regarding Item 6, the weighting w(x,y) is given by the fraction of electrons emitted from an element with the total number of emitted electrons or

wij=J(xi,yj)A(xi,yj)ijJ(xi,yj)A(xi,yj)JijAijijJijAij
(40)

if the current density does not change appreciably over the time Δt. The variation in w(x,y) was shown in Fig. 19 when J(x,y) is specified by the Fowler-DuBridge relation of Eq. (35) [which is the βF1 approximation of the right hand side of Eq. (A16)]. The weighting factor can be expressed as a “number of particles emitted from the Aij patch” nij by

nijnminwij/wmin,
(41)

where rounds to the smaller integer and nmin is the number of electrons from the patch producing the least emission. The total number of emitted electrons Ntot and average velocity components vk2 (from which “mean transverse energy” factors are evaluated) is then

Ntot=ijnij,
(42)
vk2=1Ntotijnij[vij(ta)k^]2,
(43)

where k^ is a cartesian direction vector, whereas the ij-indices refer to the grid element location.

The approximate nature of the methods meets the needs of particle-in-cell codes for an emission model that gives number, velocity, and energy spread with the rapidity of a function call. The bulk of computational effort in beam simulation codes is toward propagating the pulses through structures, and, therefore, a computationally demanding emission model is shunned. A simple model such as the Ballistic-Impulse model, if it captures the majority of the effects of surface roughness, is, therefore, desirable even at the price of limitations on its positional accuracy at the microscale or the approximate nature of the evaluation of C.

The surface angles in vk^ defined in Eq. (43) are obtained as follows. Define zs(xi,yj)zij on a grid such that xi=(i1)Δ, yj=(j1)Δ, and Δ=L/(Np1) and Np the number of grid points along an axis. Introduce two diagonal vectors,

lαΔx^+Δy^+(zi+1,j+1zi,j)z^(x^+y^+zαz^)Δ,lβΔx^Δy^+(zi,j+1zi+1,j)z^(x^y^+zβz^)Δ.
(44)

The surface normal vector n^ for Aij is then

n^ij=lα×lβlαlβ=(zα+zβ)x^+(zαzβ)y^2z^(2+zα2)(2+zβ2)
(45)

[not to be confused with the weighting factor of nij of Eq. (42)] from which the angles θij and ϕij are defined by cosθijn^ijz^ and cosϕij=(n^ijx^)/sinθ, the magnitudes of which are

cos2θij=4(2+zα2)(2+zβ2),cos2ϕij=(zα+zβ)22zα2+2zβ2+(zαzβ)2.
(46)

Alternately, if (1/2)mv(ta)2 is the kinetic energy at the anode plane, then the x^-component of the va=v(ta) is van^x and similarly for the other components. An indication of their behavior is shown in Fig. 21 using the photoemission-based choices of vo=0.3018nm/fs and κ=72.40 (parameters characteristic of photoemission as given in Table III of Ref. 21).

FIG. 21.

Trajectories (white lines) from the surface at uniformly spaced origins zs(xi,yj) evaluated until z(t)>3zc in order to make the paths clear. The background color pattern corresponds to the surface map of Fig. 16.

FIG. 21.

Trajectories (white lines) from the surface at uniformly spaced origins zs(xi,yj) evaluated until z(t)>3zc in order to make the paths clear. The background color pattern corresponds to the surface map of Fig. 16.

Close modal

Histograms of the transverse velocity component vρ2vx2+vy2 are shown in Fig. 22, where the number of times the trajectory contributed was weighted according to nij. A comparison with the similarly weighted Eq. (39) values is also shown, demonstrating the applicability of the analytic model. Such histograms can be used to generate estimates of MTE by using Eq. (39) to find the Maxwell-Boltzmann distribution factor α such that the normalized distribution f(vρ) resembles

f(vρ)2αvρexp(αvρ2),
(47)

where α may be crudely estimated by the location vρ=vmax where f(vmax) is the maximum value of the numerically evaluated ballistic-impulse histogram or α=1/(2vmax2)=4.47. A refinement to a minimizes the root mean square (rms) difference between the histogram data and Eq. (47) and gives arms=5.35. The rms fit appears as the red line in Fig. 22. Although f(vρ) may be used to estimate the MTE, using the trajectory weightings and velocity components directly to obtain

MTE12mvo2σ2sin2θ=FoL2κNtotijnijσij2(1cos2θij),
(48)

is preferable.

FIG. 22.

Normalized histogram of transverse velocity vρ=vx2+vy2 for vo=0.3018 (nm/fs) and κ=72.4 (photoemissionlike parameters). Sum over entries equals unity. “Bal-Imp” refers to the Ballistic-Impulse model of Eq. (39). “MB-fit” refers to Eq. (47) with α determined by least squares fit on the numerical data. α=5.35(fs/nm)2.

FIG. 22.

Normalized histogram of transverse velocity vρ=vx2+vy2 for vo=0.3018 (nm/fs) and κ=72.4 (photoemissionlike parameters). Sum over entries equals unity. “Bal-Imp” refers to the Ballistic-Impulse model of Eq. (39). “MB-fit” refers to Eq. (47) with α determined by least squares fit on the numerical data. α=5.35(fs/nm)2.

Close modal

The initial velocity of the trajectories is launched normal to the local surface and so does not contain velocity components parallel to the local surface. Such contributions are distinct and termed “field distortion” and “slope effects” by Bradley et al.20 As such, the normal-launch velocity factors contribute to the MTE in a manner different from “thermal emittance” for which MTEEt where E=Et+En or the sum of the transverse and normal energies, respectively. Usage of the Ballistic-Impulse model in beam optics codes will, therefore, result in larger intrinsic emittance factors associated with the cathode, compatible with the observation of Dowell and Schmerge72 that experimentally measured thermal emittance is approximately a factor of two greater than the theoretical estimate.

Although a calculation of the MTE at a plane above the rough surface is possible, it is apparent that space charge almost immediately begins to modify the trajectories of the emitted electrons. A proper evaluation of the intrinsic surface emittance will, therefore, estimate the trajectories and velocity distributions at a flat virtual cathode plane approximately at z=2zc, which can then be used as the initial conditions in a Particle-in-Cell (PIC) code, the preferred method to make predictions regarding emittance evolution.21 

The electrostatic PIC code MICHELLE73 is used to validate the SRM theory. The current density weighting factor w, field enhancement factor β, and transverse velocity distribution vρ are calculated self-consistently in MICHELLE and compared with theory. In each case, agreement between MICHELLE simulation and theory is achieved.

The mesh for the MICHELLE solver was generated using GMSH.74 The rough surface generated using SRM was calculated in MATLAB and exported on a regular grid. The regular grid surface data were mapped into GMSH by dividing each regular grid square into two triangles to create a triangulated surface. GMSH’s capability of defining a compound surface was used to combine the triangulated surface into a single geometric surface feature. In this way, an unstructured mesh of any desired mesh size can be generated in GMSH. Its application to the present surface is shown in Fig. 23.

FIG. 23.

Surface roughness model mesh vertices generated by GMSH for MICHELLE (dots) overlaid on an interpolated surface to show relative height. The surface comes from that produced by SRM and represents the rough conducting surface of the emitting region. Compare Fig. 16.

FIG. 23.

Surface roughness model mesh vertices generated by GMSH for MICHELLE (dots) overlaid on an interpolated surface to show relative height. The surface comes from that produced by SRM and represents the rough conducting surface of the emitting region. Compare Fig. 16.

Close modal

MICHELLE tracks the particle trajectories with great accuracy. In simulation, the width and length of the rough surface were normalized to one unit as was done in SRM. In simulation, an anode must exist to support the desired electrical field. The anode was positioned at a normalized distance of two units from the rough surface where the minimum z location of the rough surface is defined as zero. The four sides of the simulated vacuum space model, between the cathode and anode, are defined to be periodic such that the domain is effectively infinite in the transverse direction, and particles escaping a transverse boundary re-enter the simulation on an opposing face. Figure 24 shows a subset of the electron trajectories (green) as they emit off the rough surface and are transported toward the anode.

FIG. 24.

Particle trajectories (green) as simulated in MICHELLE. 1 in 20 particles are shown. (top) (1) Periodic boundary. (2) Cathode with rough surface. (3) Vacuum space and trajectories. (4) Vacuum space outline (cut at y=0). (bottom) Cathode surface showing grid.

FIG. 24.

Particle trajectories (green) as simulated in MICHELLE. 1 in 20 particles are shown. (top) (1) Periodic boundary. (2) Cathode with rough surface. (3) Vacuum space and trajectories. (4) Vacuum space outline (cut at y=0). (bottom) Cathode surface showing grid.

Close modal

The field enhancement factor β is evaluated by the surface field normalized against the bulk background field. Figure 25 shows the results of the MICHELLE simulation and its agreement with SRM. The peak enhancement in SRM was approximately 1.5 and in MICHELLE, 1.41. The current density weighting factor w is normalized against the minimum emitted current density. For the given rough surface, the minimum current density occurs near (0, 0) where there is a large depression in the surface shielding the cathode from the electric field, resulting in reduced emission. The current density weighting factor as simulated in MICHELLE is shown in Fig. 26. To first order, SRM and MICHELLE agree. The peak weighting factor in SRM is higher. This is due to the normalization against the minimum field which is very sensitive to the deep depression in the middle of the rough surface. Careful meshing of this region is needed to properly resolve the low field in the surface depression.

FIG. 25.

Field enhancement factor β(x,y,zs) using MICHELLE on the normalized geometry is evaluated at the black dots; the magnitude of β at them is conveyed by the interpolated color map. Compare Fig. 17.

FIG. 25.

Field enhancement factor β(x,y,zs) using MICHELLE on the normalized geometry is evaluated at the black dots; the magnitude of β at them is conveyed by the interpolated color map. Compare Fig. 17.

Close modal
FIG. 26.

Weighting factor w(x,y,zs) for J/JFD [see Eq. (35)] using MICHELLE is evaluated at the black dots; the magnitude of w at them is conveyed by the interpolated color map. Compare Fig. 19. Parameters: Φ=2.3eV, T=300K, and λ=532nm. The number of electrons emitted from each area element A(x;y) is proportional to w(x;y).

FIG. 26.

Weighting factor w(x,y,zs) for J/JFD [see Eq. (35)] using MICHELLE is evaluated at the black dots; the magnitude of w at them is conveyed by the interpolated color map. Compare Fig. 19. Parameters: Φ=2.3eV, T=300K, and λ=532nm. The number of electrons emitted from each area element A(x;y) is proportional to w(x;y).

Close modal

Finally, the transverse velocity of the particles was modeled in MICHELLE and compared with SRM. In SRM, the average electron energy at the z location of the calculation was 64 eV. The anode potential in MICHELLE was adjusted such that the average potential at the export location was 64 V. The simulated velocity distribution is in agreement with SRM as can be seen by comparing the qualitative behavior shown in Figs. 22–27.

FIG. 27.

Normalized histogram of transverse velocity vρ=vx2+vy2 for vo=0.3018nm/fs as simulated in MICHELLE. Compare to Fig. 22.

FIG. 27.

Normalized histogram of transverse velocity vρ=vx2+vy2 for vo=0.3018nm/fs as simulated in MICHELLE. Compare to Fig. 22.

Close modal

A theory and methodology able to represent surface features, from an ordered array of emitters to the random undulations characteristic of cathode surface roughness, was developed by combining a general thermal-field-photoemission theory with the point charge model and ballistic-impulse methods for launching electrons. Field enhancement effects and screening (aka “shielding”) depend on a range of feature size scales and affect what launch velocities are appropriate in determining particle trajectories from surfaces. The methodology enables accounting for structure at the scale of emission sites (“microscale”) which is generally difficult to resolve, up to feature sizes characteristic of the pitch of an ordered array or visible in scanning electron microscope (SEM) imaging of surface features (“mesoscale”). As the field variation along surface features as well as along electron trajectory lines can be resolved, simultaneous predictions of thermal-field emission contributions, local current density or photoemission yield, and velocity distributions can be made.

The methodology was (i) applied to square and triangular arrays of conical emitters characteristic of field and thermal-field sources showing the onset of tunneling effects in thermal-field emission and its dependence on pitch and arrangement; (ii) applied to a surface roughness model showing the ability to handle random variation effects; (iii) used to show how work function characterization is affected by surface structure for thermal emitters; (iv) applied to show how conventional field enhancement, work function, and emission area estimates are undermined; (v) used to show how unit cell models can be used to give velocity distributions and trajectories needed by PIC codes that mitigate difficult grid size resolution demands related to microscale variation; and (vi) used to show how semianalytic predictions of MTE factors may be approximated to give a measure of intrinsic emittance at the cathode due to extended arrays of emitters or random surface roughness effects.

Two configurations were considered in depth: a periodic array with nominally identical emitters and a surface with random variations in geometry. With regard to the former, the methodology was applied to a geometry fashioned after a conical (C12A7) emitter: it was shown that thermal-field effects can have a pronounced impact on the emission parameters that are ascribed to emitters within an array based on current-voltage data, wherein the putative values depart substantially from the actual ones of the theoretical model. With regard to the latter, a revised PCM was applied to surface roughness, and random variation in the surface was shown to contribute to transverse velocity components at the rough surface: it was shown that an updated ballistic-impulse model performs well in predicting the velocity components at a “virtual anode” plane above the surface of the cathode at a distance approximately twice the surface roughness scale. The two cases have, therefore, vetted the methodology as a means to include microscale emission effects in a mesoscale framework suitable for beam optics codes. To confirm that, simulations comparing the theoretical models as contained in the code MICHELLE demonstrated close agreement with the analytical and seminumerical models applied to surface roughness.

The authors gratefully acknowledge support from the following sources: (i) the Air Force Office of Scientific Research (AFOSR) Plasma and Electro-Energetic Physics Program under Contract No. 18RDCOR016 and (ii) the Laboratory Directed Research and Development (LDRD) program of Los Alamos National Laboratory (LANL) under Project No. 20150394DR (LANL is operated by Los Alamos National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy, Contract No. DE-AC52-06NA25396). Additionally, the authors thank J. L. Shaw (NRL) for discussions regarding characterizations of C12A7 emitter cones.

An algorithm to find the factors n(F,T) and s(F,T) needed by Eq. (16) for the evaluation of the general thermal-field and photoemission (GTFP) current density is given, synopsizing prior treatments22–24 for numerical efficiency. The compact relation to the canonical FN, RLD, and FD (“canonical”) equations is as in Appendix A of Ref. 21.

The thermal (βT) and field (βF) slope factors are obtained from the gradients with energy of the exponents of the Fermi-Dirac distribution (thermal) and the Gamow factor θ(E) on which transmission probability is based (field). The energy slope factor βT for temperature is independent of energy and given by

βT=E[(μE)/kBT]=1/kBT,
(A1)

where μ is the Fermi energy. The energy slope factor βF(E) for field depends on energy and is related to the gradient of the Gamow factor θ(E) by

βF(Em)Eθ(E)|E=Em,
(A2)

where Em is the location of the maximum of the product of the transmission probability and the supply function [D(E)f(E)], and the Gamow factor is defined by

θ(E)22mxx+U(x)Edx,
(A3)

where U(x) is the one-dimensional potential energy including the image charge modification, U(x±)E=0, and E in the one-dimensional formulation corresponds to the normal energy in a three-dimensional. The resulting potential energy barrier, aka the Schottky-Nordheim (SN) barrier, subject to an electric field E, is

U(x)=μ+Φq|E|xq216πε0xμ+ΦFxQx,
(A4)

where F and Q are defined by these equations. The Gamow factor θ(E) is taken as a cubic in energy for energies above the Fermi energy and below the barrier maximum,22 that is, for when μEμ+Φ4QF, which implies that βF(E) is a quadratic in energy. The coefficients of the interpolating polynomial can be inferred from the behavior of θ(E) near E=μ (which is well approximated by the Fowler-Nordheim representation) and near E=μ+Φ4QF (which is well approximated by the parabolic barrier representation). The Fowler-Nordheim representation introduces two functions v(y) and t(y), sometimes called “Schottky-Nordheim barrier functions,” which depend on elliptical integrals26 but for which very good approximations are75 

v(y)=1y23(3ln(y)),
(A5)
t(y)=1+y29(1ln(y)),
(A6)

where y4QF/Φ and is related to the Schottky-reduced barrier term,

ϕΦ4QF=(1y)Φ.
(A7)

Coefficients of the cubic polynomial can be uniquely determined by specifying the values θ(μ) and θ(μ+ϕ) and derivatives βF(μ) and βF(μ+ϕ) at the end points, which are known to be62,76

θF(μ)=43F2mΦ3v(y),
(A8)
βF(μ)=2F2mΦt(y),
(A9)
θF(μ+ϕ)=0,
(A10)
βF(μ+ϕ)=πFmΦy.
(A11)

An expedient algorithm to find JGTF(F,T) is then given by introducing the temperatures

Tmin=1kBβF(μ),Tmax=1kBβF(μ+ϕ)
(A12)

and then evaluating n(F,T) and s(F,T) according to

  • For T<Tmin:

    set n=βT/βF(μ) and s=θ(μ)

  • For T>Tmax:

    set n=βT/βF(μ+ϕ) and s=βF(μ+ϕ)ϕ

  • Otherwise:

    set n=1, let Em be defined by βT=βF(Em) and
    s=βT(Em+θ(Em)βF(Em)μ).
    (A13)

Construct

JGTFP(F,T)ARLDT2N(n,s),
(A14)

where ARLD=(qmkB2)/(2π23)=120.1735A/cm2K2. Let N(n,s) be defined by

N(n,s)=nln[1+en(zs)]1+ezdz.
(A15)

It can be shown that

N(n,s)+N(n,s)=12n2s2+ζ(2)(n2+1),
(A16)
N(1,s)=(s+1)es,
(A17)

where ζ(2)=π2/6 (Riemann Zeta function). A useful approximation is enabled by introducing

Σ(x)1+x21x2+(π262)x2+(7π43602)x4
(A18)

in terms of which

N(n,s)esn2Σ(1n)+ensΣ(n).
(A19)

The following limits represent the connection with the canonical equations:23 

limF0JGTF(F,T)=ARLDT2exp(ϕ/kBT),
(A20)
limT0JGTF(F,T)=ARLD(kBβF(μ))2exp(θ(μ)),
(A21)

which is tantamount to n0 and n, respectively. Equation (A20) is the Richardson-Laue-Dushman (RLD) equation25 when small reflection coefficients are neglected and an “effective” work function ϕ=Φ4QF is used,77 and Eq. (A21) is the Murphy and Good26 form of the Fowler-Nordheim equation when Eqs. (A8) and (A9) are used. Cases are shown in Figs. 28 and 29. Particular numerical cases are in Table III, in which conditions are shown where the RLD and FN equations both underpredict emitted current, thereby having an impact on predicted Φ and β values derived from I(V) data.

FIG. 28.

JGTF(F,T) as a function of F for T=1173K and Φ=4.5eV. Also shown are Eq. (A20) (red) and Eq. (A21) (blue) and gray lines corresponding to F(Tmin)=1.361eV/nm and F(Tmax)=2.273eV/nm as per Eq. (A12).

FIG. 28.

JGTF(F,T) as a function of F for T=1173K and Φ=4.5eV. Also shown are Eq. (A20) (red) and Eq. (A21) (blue) and gray lines corresponding to F(Tmin)=1.361eV/nm and F(Tmax)=2.273eV/nm as per Eq. (A12).

Close modal
FIG. 29.

JGTF(F,T) as a function of F for T=1173K and Φ=2.1eV. Also shown are Eq. (A20) (red) and Eq. (A21) (blue) and gray lines corresponding to F(Tmin)=1.36eV/nm and F(Tmax)=1.617eV/nm as per Eq. (A12).

FIG. 29.

JGTF(F,T) as a function of F for T=1173K and Φ=2.1eV. Also shown are Eq. (A20) (red) and Eq. (A21) (blue) and gray lines corresponding to F(Tmin)=1.36eV/nm and F(Tmax)=1.617eV/nm as per Eq. (A12).

Close modal
TABLE III.

Evaluations of n, s, and J for various fields F (eV/nm), temperature T (K), and work functions Φ (eV). Thermal-field current densities J are in A/cm2.

FTΦnslogJGTFlogJRLDlogJFN
0.01 1400 2.0 0.021 741.173 1.605 1.604  
2.00 1100 4.5 1.000 26.265 −1.809 −4.680 −3.381 
2.00 1100 2.0 1.332 2.324 7.738 6.775 7.402 
2.00 1400 4.5 1.000 22.993 −0.234 −1.718 −3.381 
2.00 1400 2.0 1.047 2.324 7.899 7.281 7.402 
5.00 300 4.5 8.395 7.605 5.589 −23.486 5.579 
5.00 1400 4.5 1.799 7.605 5.827 1.832 5.579 
8.00 77 4.5 51.311 3.078 7.937  7.936 
FTΦnslogJGTFlogJRLDlogJFN
0.01 1400 2.0 0.021 741.173 1.605 1.604  
2.00 1100 4.5 1.000 26.265 −1.809 −4.680 −3.381 
2.00 1100 2.0 1.332 2.324 7.738 6.775 7.402 
2.00 1400 4.5 1.000 22.993 −0.234 −1.718 −3.381 
2.00 1400 2.0 1.047 2.324 7.899 7.281 7.402 
5.00 300 4.5 8.395 7.605 5.589 −23.486 5.579 
5.00 1400 4.5 1.799 7.605 5.827 1.832 5.579 
8.00 77 4.5 51.311 3.078 7.937  7.936 

For photoemission, the factor (μ+ϕE) appearing in the transmission probability has E augmented by the photon energy ω which renders (μ+ϕEω)<0. Under typical photoemission conditions, N(n,s) now becomes negligible compared to N(n,s), and so the emitted current density is proportional to23,24

T2N(n,s)12(ωϕ)2+ζ(2)(βT2+βF2).
(A22)

When fields are small, then βFβT and the temperature-dependent Fowler-DuBridge equation78 [see Eq. (35)] results.

The point charge model gives an analytic means to determine the field enhancement and radius at the apex of an isolated emitter characterized by two factors n (number of point charges) and r (governing how the height of the emitter increases with each point charge).57 Introduce the matrix M and vectors S and P for the dipole defined by55 

[M]nj2Sj(r)Sn+1j(r)[Sj(r)+Sn+1(r)]θnj,[S]jSj+1(r),[P]jPj(r),
\vskip 7.6pt(B1)

where θnj=1 for jn and 0 otherwise, and the prime on S reinforces the j+1 index in its definition. The vector P is the solution to the matrix equation,

MP=S.
(B2)

The element Pj+1 is directly expressible in terms of the lower index Pj, possible because M is a lower-triangular matrix. Similarly, let Qj=Pj/(rjSj). Then, the βn is found from

βn=1+j=1nSn+1QjMnj2.
(B3)

Let Tj=Qj/(4rjSj). Then, the an is found from

an=j=1nSn+1QjMnj2j=1nTj(Sj2+3Sn+12)Mnj3.
(B4)

Observe that a factor of Mnj3 is present in the denominator of Eq. (B4): this corrects Eq. (15) of Ref. 55, which mistakenly represented the power as Mnj2. Although not presented as such here, Eqs. (B3) and (B4) may be cast as matrix and vector multiplications.

Example values are given in Table IV for the particular case r=0.7. Observe that to a reasonable approximation, βn/βN1.24nN, reflecting Schottky’s conjecture regarding the multiplicity of β-factors.37 Similarly, aN/an0.72Nn for the chosen value of r. The exponential behavior of both is shown in Fig. 30.

FIG. 30.

Demonstration of exponential behavior, implying Schottky’s conjecture, of apex beta factor βn and apex radius an for N=10 and r=0.7.

FIG. 30.

Demonstration of exponential behavior, implying Schottky’s conjecture, of apex beta factor βn and apex radius an for N=10 and r=0.7.

Close modal
TABLE IV.

Values of Sn(r) [Eq. (2)], Pn(r) [Eq. (B2)], βn(r) [Eq. (B3)], and an(r) [Eq. (B4)] evaluated for r=0.7.

nSn(r)Pn(r)βn(r)βn/βn1an(r)
1.0000 2.295 0 4.0582  0.664 53 
1.7000 1.537 2 5.0640 1.2479 0.557 21 
2.1900 1.308 5 6.3931 1.2625 0.420 66 
2.5330 1.131 2 8.0182 1.2542 0.308 62 
2.7731 0.979 60 9.9999 1.2472 0.223 29 
2.9412 0.847 85 12.424 1.2424 0.160 19 
3.0588 0.733 29 15.398 1.2393 0.114 27 
3.1412 0.633 87 19.053 1.2374 0.081 171 
3.1988 0.547 74 23.553 1.2362 0.057 481 
10 3.2392 0.473 22 29.098 1.2354 0.040 608 
nSn(r)Pn(r)βn(r)βn/βn1an(r)
1.0000 2.295 0 4.0582  0.664 53 
1.7000 1.537 2 5.0640 1.2479 0.557 21 
2.1900 1.308 5 6.3931 1.2625 0.420 66 
2.5330 1.131 2 8.0182 1.2542 0.308 62 
2.7731 0.979 60 9.9999 1.2472 0.223 29 
2.9412 0.847 85 12.424 1.2424 0.160 19 
3.0588 0.733 29 15.398 1.2393 0.114 27 
3.1412 0.633 87 19.053 1.2374 0.081 171 
3.1988 0.547 74 23.553 1.2362 0.057 481 
10 3.2392 0.473 22 29.098 1.2354 0.040 608 

In the spherical coordinate system where (ρ/a)2+(z/a)2=1, then

ρ=asinθ,z=acosθ.
(C1)

For the prolate spheroidal coordinate system,62 then

ρ=Lsinhηsinυ,z=Lcoshηcosυ.
(C2)

The ellipsoidal coordinates are defined by (z/be)2+(ρ/ae)2=1 with ae=Lsinhη and be=Lcoshη in Eq. (C2) so that

ρ=aesinυ,z=becosυ.
(C3)

The spherical case is recovered in the limit υ for which ae/be1.

The hyperbolic coordinates are defined by (z/bh)2(ρ/ah)2=1 with ah=Lsinυ and bh=Lcosυ in Eq. (C2) so that

ρ=ahsinhη,z=bhcoshη.
(C4)

The rewriting of Eqs. (C1), (C3), and (C4) along the x-axis (for which ρ=x) for two pairs of points (x1,z1) at the apex and (x2,z2) off axis gives rise to Eq. (12) after correcting for the placement of the origin on the z-axis.

The impulse model, introduced in Ref. 21, has been substantially revised to enable a fast analytical algorithm designed for usage in beam optics codes such as MICHELLE.73,79 The model is most easily demonstrated by considering trajectories from a basic hemispherical protrusion. The equations are dimensioned for position and velocity, and so the length scale a (radius of the hemisphere) and factor Fo (product of unit charge q and background electric field magnitude) are explicitly shown. The potential energy associated with metal hemisphere is defined in spherical coordinates (r,θ) with boundary conditions such U(a,θ)=0, giving

U(r,θ)=Forcosθ[1(ar)3].
(D1)

The image charge contribution is ignored, being significant for only nanometers from the surface. The radial gradient of the potential energy gives force as

F(r,θ)=rU(r,θ)=Focosθ[1+2(ar)3],
(D2)

where Fo is taken to be the asymptotic (or far) field: when evaluated at the surface (r=a), the field enhancement factor defined by F(a,θ)=β(θ)Fo is then seen to be β(θ)=3cosθ and is independent of the radius, although that is not true in general for rough surfaces. The field components in cylindrical coordinates defined by r2=ρ2+z2 and tanθ=ρ/z (where z and ρ have dimensions of length while finding the ballistic equations) are given by

Fρ(ρ,z)=3Foa3ρzr5,Fz(ρ,z)=3Foa3z2r5+Fo[1(ar)3].
(D3)

The enhanced field near the boss operates only for a short time before the trajectory is governed by the constant field Fo, e.g., for ρ=0 and z/a=(2,3), Fz/Fo=(1.25,1.07).

Letting ρ2=x2+y2, then from finite difference methods, the following formulae allow for finding the trajectories ρ(t) and z(t) by

ρ(t+δt)ρ(t)+vρ(t)δt+Fρ[ρ(t),z(t)]2mδt2=ρ(t)+12[vρ(t+δt)+vρ(t)]δt,z(t+δt)z(t)+vz(t)δt+Fz[ρ(t),z(t)]2mδt2
(D4)
=z(t)+12[vz(t+δt)+vz(t)]δt.
(D5)

These relations are exact if the forces Fρ and Fz are constant: the resulting equations will then be referred to as “ballistic.” Under ballistic conditions,

ρ(t)=ρ(0)+12[vρ(t)+vρ(0)]t,
(D6)
z(t)=z(0)+12[vz(t)+vz(0)]t.
(D7)

The ballistic equations are then tantamount to Fρ0 and FzFo, with the launch sites defined by z(0)=acosθ and ρ(0)=asinθ. Introduce pvot/a as a dimensionless time and κ as the ratio of the energy scale associated with the background field and that associated with the initial kinetic energy defined by

κFoamvo2=Foτmvo,
(D8)

where τ=a/vo, then the ballistic trajectories are defined by (ρb(t),zb(t)), and are given by

ρb(t)a=(1+p)sinθ,
(D9)
zb(t)a=(1+p)cosθ+12κp2.
(D10)

The “Impulse Approximation” augments the initial velocity voσvo to account for the short duration when the accelerating field is not constant and thereby regains usage of the ballistic relations away from roughness effects. The revised impulse theory now does so by making the replacement

vρ(t)vρ(0)+1mFρ[ρ(0),z(0)]t,
(D11)
vz(t)vz(0)+1mFz[ρ(0),z(0)]t.
(D12)

Using Eqs. (D3) and (D11) in Eq. (D6) results in

ρ(t)a=sinθ{1+p[1+32κpcosθ]},
(D13)
z(t)a=cosθ{1+p[1+32κpcosθ]}+12κp2.
(D14)

Assuming the impulse occurs only close to the hemisphere for a time p=a/vot1 identifies the impulse factor σ as

σ(θ)1+32Cκcosθ,
(D15)

where C is a constant of order unity chosen to give the best agreement between the impulse trajectories and the actual trajectories. It may be crudely estimated by demanding that on axis (θ=0), a particle passing the location z=na for n2 has the same energy kinetic (1/2)mv(t)2 for both the Ballistic-Impulse approximation and the actual trajectory solution, that is, from the equations relating changes of kinetic energy to changes in potential energy

12mv212mvo2=U(0,a)U(0,na),
(D16)
12mv212m(σvo)2=Foa(Fona).
(D17)

Substituting the first into the second equation, letting σ1+(3/2)Cκ (on axis), assuming n21, and collecting gives

C43(1+1+2κ).
(D18)

Such a relationship, although approximately adequate, should not take precedence over C determined numerically, as in Fig. 31, where C is adjusted to provide good correspondence: two examples are shown with symbols referring to numerical evaluation of the trajectories by Eq. (D15), and lines for using Eqs. (D13) and (D14) for the value of C shown [which are somewhat higher than the predictions of Eq. (D18)]. The revised Ballistic-Impulse model performs better than prior theory (compare Figs. 12–14 of Ref. 21) and is better suited to approximating the effects of surface roughness.

FIG. 31.

Numerically calculated trajectories (symbol) compared to Ballistic-Impulse trajectories (red lines) for the values of κ and C shown; note that Eq. (D18) would give C(κ=1)=0.488 and C(κ=10)=0.239. The new Ballistic-Impulse model is also seen to provide a better account of the trajectories than the approximations of Ref. 21.

FIG. 31.

Numerically calculated trajectories (symbol) compared to Ballistic-Impulse trajectories (red lines) for the values of κ and C shown; note that Eq. (D18) would give C(κ=1)=0.488 and C(κ=10)=0.239. The new Ballistic-Impulse model is also seen to provide a better account of the trajectories than the approximations of Ref. 21.

Close modal
1.
J. R.
Harris
,
K. L.
Jensen
,
J. J.
Petillo
,
S.
Maestas
,
W.
Tang
, and
D. A.
Shiffler
,
J. Appl. Phys.
121
,
203303
(
2017
).
2.
R. U.
Martinelli
,
Appl. Opt.
12
,
1841
(
1973
).
3.
C.
Ducati
,
E.
Barborini
,
P.
Piseri
,
P.
Milani
, and
J.
Robertson
,
J. Appl. Phys.
92
,
5482
(
2002
).
4.
F. M.
Charbonnier
,
C.
Bennette
, and
L.
Swanson
,
J. Appl. Phys.
38
,
627
(
1967
).
5.
D.
Wright
and
J.
Woods
,
Proc. Phys. Soc. Lond. B
65
,
134
(
1952
).
6.
J.
Robertson
, “Theory of electron field emission from diamond and diamond-like carbon,” in Materials Research Society Symposium Proceedings (Materials Research Society, Warrendale, PA, 1998), Vol. 498, pp. 197–208.
7.
D.
Chen
,
R.
Jacobs
,
V.
Vlahos
,
K. L.
Jensen
,
D.
Morgan
, and
J.
Booske
, in IEEE International Vacuum Electronics Conference (IVEC) (IEEE, 2018), pp. 39–40.
8.
H. C.
Liu
and
D. D.
Coon
,
J. Appl. Phys.
64
,
6785
(
1988
).
9.
J. M.
Finn
,
T. M.
Antonsen
, and
W. M.
Manheimer
,
IEEE Trans. Plasma Sci.
16
,
281
(
1988
).
10.
Y. Y.
Lau
,
J. Appl. Phys.
61
,
36
(
1987
).
11.
J.
Zuber
,
K. L.
Jensen
, and
T.
Sullivan
,
J. Appl. Phys.
91
,
9379
(
2002
).
12.
R.
Miller
,
Y. Y.
Lau
, and
J. H.
Booske
,
Appl. Phys. Lett.
91
,
074105
(
2007
).
13.
R.
Miller
,
Y.
Lau
, and
J.
Booske
,
J. Appl. Phys.
106
,
104903
(
2009
).
14.
J. R.
Harris
,
D. A.
Shiffler
,
K. L.
Jensen
, and
J. W.
Lewellen
,
J. Appl. Phys.
125
,
215307
(
2019
).
15.
M.
Krasilnikov
, “THPPH013: Impact of the cathode roughness on the emittance of an electron beam,” in Proceedings of FEL Conference (BESSY, Berlin, Germany, 2006), pp. 583–586.
16.
S.
Karkare
and
I. V.
Bazarov
,
Appl. Phys. Lett.
98
,
094104
(
2011
).
17.
T.
Vecchione
,
J.
Feng
,
W.
Wan
,
H. A.
Padmore
,
I.
Ben-Zvi
,
X.
Liang
,
M.
Ruiz-Oses
,
T.
Rao
,
J.
Smedley
, and
D.
Dowell
, in International Particle Accelerator Conference (IEEE, New Orleans, LA, 2012), p. MOPPP041.
18.
T.
Charles
,
D.
Paganin
, and
R.
Dowd
,
Nucl. Instrum. Methods Phys. Res. A
828
,
201
(
2016
).
19.
D. A.
Dimitrov
,
G. I.
Bell
,
J.
Smedley
,
I.
Ben-Zvi
,
J.
Feng
,
S.
Karkare
, and
H. A.
Padmore
,
J. Appl. Phys.
122
,
165303
(
2017
).
20.
D. J.
Bradley
,
M. B.
Allenson
, and
B. R.
Holeman
,
J. Phys. D Appl. Phys.
10
,
111
(
1977
).
21.
K. L.
Jensen
,
D. A.
Shiffler
,
J. J.
Petillo
,
Z.
Pan
, and
J. W.
Luginsland
,
Phys. Rev. Spec. Top. Accel. Beams
17
,
043402
(
2014
).
22.
K. L.
Jensen
,
J. Appl. Phys.
102
,
024911
(
2007
).
23.
K. L.
Jensen
,
Introduction to the Physics of Electron Emission
(
John Wiley & Sons, Inc.
,
Hoboken, NJ
,
2017
).
24.
K. L.
Jensen
,
IEEE Trans. Plasma Sci.
46
,
1881
(
2018
).
25.
C.
Herring
and
M.
Nichols
,
Rev. Mod. Phys.
21
,
185
(
1949
).
26.
E. L.
Murphy
and
R. H.
Good
,
Phys. Rev.
102
,
1464
(
1956
).
27.
R. A.
Millikan
and
C.
Lauritsen
,
Proc. Natl. Acad. Sci. U.S.A.
14
,
45
(
1928
).
28.
W. W.
Dolan
and
W. P.
Dyke
,
Phys. Rev.
95
,
327
(
1954
).
29.
J. R.
Harris
,
K. L.
Jensen
, and
D. A.
Shiffler
,
J. Phys. D Appl. Phys.
48
,
385203
(
2015
).
30.
L. D.
Filip
,
J. D.
Carey
, and
S. R. P.
Silva
,
J. Appl. Phys.
109
,
084527
(
2011
).
31.
L.
Eyges
,
The Classical Electromagnetic Field
(
Addison-Wesley Publishing Company
,
Reading, MA
,
1972
), p. 62.
32.
J.
He
,
P.
Cutler
,
N.
Miskovsky
,
T.
Feuchtwang
,
T.
Sullivan
, and
M.
Chung
,
Surf. Sci.
246
,
348
(
1991
).
33.
D.
Kirkpatrick
,
A.
Mankofsky
, and
K.
Tsang
,
Appl. Phys. Lett.
60
,
2065
(
1992
).
34.
A.
Kusne
and
D.
Lambeth
,
IEEE Trans. Electron Devices
57
,
712
(
2010
).
35.
M.
Pant
and
L. K.
Ang
,
J. Appl. Phys.
107
,
106103
(
2010
).
36.
T.
Stern
,
B.
Gossling
, and
R.
Fowler
,
Proc. R. Soc. Lond. A
124
,
699
(
1929
).
37.
K. L.
Jensen
,
D. A.
Shiffler
,
J. R.
Harris
, and
J. J.
Petillo
,
AIP Adv.
6
,
065005
(
2016
).
38.
T.
Tsang
,
T.
Srinivasanrao
, and
J.
Fischer
,
Nucl. Instrum. Methods A
318
,
270
(
1992
).
39.
L.
Ang
,
Y.
Lau
,
R.
Gilgenbach
, and
H.
Spindler
,
Appl. Phys. Lett.
70
,
696
(
1997
).
40.
J.
Smedley
,
T.
Rao
, and
Q.
Zhao
,
J. Appl. Phys.
98
,
043111
(
2005
).
41.
D.
Dowell
,
I.
Bazarov
,
B.
Dunham
,
K.
Harkay
,
C.
Hernandez-Garcia
,
R.
Legg
,
H.
Padmore
,
T.
Rao
,
J.
Smedley
, and
W.
Wan
,
Nucl. Instrum. Methods Phys. Res. A
622
,
685
(
2010
).
42.
D. A.
Dimitrov
,
D.
Smithe
,
J. R.
Cary
,
I.
Ben-Zvi
,
T.
Rao
,
J.
Smedley
, and
E.
Wang
,
J. Appl. Phys.
117
,
055708
(
2015
).
43.
N. A.
Moody
,
K. L.
Jensen
,
A.
Shabaev
,
S. G.
Lambrakos
,
J.
Smedley
,
D.
Finkenstadt
,
J. M.
Pietryga
,
P. M.
Anisimov
,
V.
Pavlenko
,
E. R.
Batista
,
J. W.
Lewellen
,
F.
Liu
,
G.
Gupta
,
A.
Mohite
,
H.
Yamaguchi
,
M. A.
Hoffbauer
, and
I.
Robel
,
Phys. Rev. Appl.
10
,
047002
(
2018
).
44.
D. H.
Dowell
,
F. K.
King
,
R. E.
Kirby
,
J. F.
Schmerge
, and
J. M.
Smedley
,
Phys. Rev. Spec. Top. Accel. Beams
9
,
063502
(
2006
).
45.
S.
Matsuishi
,
Y.
Toda
,
M.
Miyakawa
,
K.
Hayashi
,
T.
Kamiya
,
M.
Hirano
,
I.
Tanaka
, and
H.
Hosono
,
Science
301
,
626
(
2003
).
46.
Y.
Toda
,
H.
Yanagi
,
E.
Ikenaga
,
J.
Kim
,
M.
Kobata
,
S.
Ueda
,
T.
Kamiya
,
M.
Hirano
,
K.
Kobayashi
, and
H.
Hosono
,
Adv. Mater.
19
,
3564
(
2007
).
47.
L. P.
Rand
and
J. D.
Williams
,
IEEE Trans. Plasma Sci.
43
,
190
(
2015
).
48.
E.
Feizi
and
A. K.
Ray
,
J. Display Technol.
12
,
451
(
2016
).
49.
F.
Jin
,
Y.
Liu
, and
C. M.
Day
,
Appl. Phys. Lett.
88
,
163116
(
2006
).
50.
V. V.
Zhirnov
,
E. I.
Givargizov
, and
P. S.
Plekhanov
,
J. Vac. Sci. Technol. B
13
,
418
(
1995
).
51.
Y.
Toda
,
S.
Matsuishi
,
K.
Hayashi
,
K.
Ueda
,
T.
Kamiya
,
M.
Hirano
, and
H.
Hosono
,
Adv. Mater.
16
,
685
(
2004
).
52.
Y.
Toda
,
S. W.
Kim
,
K.
Hayashi
,
M.
Hirano
,
T.
Kamiya
,
H.
Hosono
,
T.
Haraguchi
, and
H.
Yasuda
,
Appl. Phys. Lett.
87
,
254103
(
2005
).
53.
S.
Jian
,
Z.
Dan
,
Z.
Fei-Hu
, and
G.
Yang
,
Appl. Surf. Sci.
433
,
358
(
2018
).
54.
K. L.
Jensen
,
Y. Y.
Lau
,
D. W.
Feldman
, and
P. G.
O’Shea
,
Phys. Rev. Spec. Top. Accel. Beams
11
,
081001
(
2008
).
55.
K. L.
Jensen
,
J. Appl. Phys.
107
,
014905
(
2010
).
56.
The introduction of the notation lo to describe the characteristic length scale is a small notational departure from prior treatments in the interest of clarity.
57.
K. L.
Jensen
,
P. G.
O’Shea
,
D. W.
Feldman
, and
J. L.
Shaw
,
J. Appl. Phys.
107
,
014903
(
2010
).
58.
J. R.
Harris
,
K. L.
Jensen
,
D. A.
Shiffler
, and
J. J.
Petillo
,
Appl. Phys. Lett.
106
,
201603
(
2015
).
59.
R. G.
Forbes
,
Appl. Phys. Lett.
110
,
133109
(
2017
).
60.
G. D.
Smith
, Numerical Solution of Partial Differential Equations: Finite Difference Methods, Oxford Applied Mathematics and Computing Science Series (Clarendon Press, Oxford University Press, Oxford, UK, New York, 1985).
61.
T. A.
de Assis
and
F. F.
Dall’Agnol
,
Nanotechnology
27
,
44LT01
(
2016
).
62.
K. L.
Jensen
,
D. A.
Shiffler
,
M.
Peckerar
,
J. R.
Harris
, and
J. J.
Petillo
,
J. Appl. Phys.
122
,
064501
(
2017
).
63.
A. C.
Keser
,
T. M.
Antonsen
,
G. S.
Nusinovich
,
D. G.
Kashyn
, and
K. L.
Jensen
,
Phys. Rev. Spec. Top. Accel. Beams
16
,
092001
(
2013
).
64.
K. L.
Jensen
, in Wiley Encyclopedia of Electrical and Electronics Engineering, edited by J. G. Webster (John Wiley & Sons, Inc., New York, 2014).
65.
An erroneous factor of 1 in the denominator in Eq. (105) of Ref. 64 has been corrected in Eq. (25).
66.
L.
Nilsson
,
O.
Groening
,
O.
Kuettel
,
P.
Groening
, and
L.
Schlapbach
,
J. Vac. Sci. Technol. B
20
,
326
(
2002
).
67.
D. A.
Shiffler
,
W.
Tang
,
K. L.
Jensen
,
K.
Golby
,
M.
LaCour
,
J. J.
Petillo
, and
J. R.
Harris
,
J. Appl. Phys.
118
,
083302
(
2015
).
68.
J. W.
Lewellen
and
C. A.
Brau
,
Nucl. Instrum. Methods Phys. Res. A
507
,
323
(
2003
).
69.
G.
Gevorkyan
,
S.
Karkare
,
S.
Emamian
,
I. V.
Bazarov
, and
H. A.
Padmore
,
Phys. Rev. Accel. Beams
21
,
093401
(
2018
).
70.
J. R.
Harris
and
K. L.
Jensen
,
J. Appl. Phys.
125
,
215306
(
2019
).
71.
Recall that because (x,y,z) are dimensionless, the factor of L that attends z~=L1z must be retained, where z~ is the dimensioned length.
72.
D. H.
Dowell
and
J. F.
Schmerge
,
Phys. Rev. Spec. Top. Accel. Beams
12
,
074201
(
2009
).
73.
J.
Petillo
,
K.
Eppley
,
D.
Panagos
,
P.
Blanchard
,
E.
Nelson
,
N.
Dionne
,
J.
DeFord
,
B.
Held
,
L.
Chernyakova
,
W.
Krueger
,
S.
Humphries
,
T.
McClure
,
A.
Mondelli
,
J.
Burdette
,
M.
Cattelino
,
R.
True
,
K.
Nguyen
, and
B.
Levush
,
IEEE Trans. Plasma Sci.
30
,
1238
(
2002
).
74.
C.
Geuzaine
and
J.-F.
Remacle
,
Int. J. Numer. Methods Eng.
79
,
1309
(
2009
).
75.
J. H. B.
Deane
and
R. G.
Forbes
,
J. Phys. A Math. Theor.
41
,
395301
(
2008
).
76.
K. L.
Jensen
,
J. Appl. Phys.
111
,
054916
(
2012
).
77.
We thank an anonymous reviewer for bringing our attention to Eq. (17) of K. T. Compton and I. Langmuir, Rev. Mod. Phys. 2, 123–242 (1930) who give the RLD equation [their Eq. (6)] in the form we have, identify ϕ=Φ4QF as the “effective” work function as we have, and use F=q|E| as we do.
78.
L. A.
DuBridge
,
Phys. Rev.
43
,
0727
(
1933
).
79.
K. L.
Jensen
,
J. J.
Petillo
,
S.
Ovtchinnikov
,
D. N.
Panagos
,
N. A.
Moody
, and
S. G.
Lambrakos
,
J. Appl. Phys.
122
,
164501
(
2017
).