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.

## I. INTRODUCTION

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 $\mu 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 breakdown^{4} 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 yield^{15–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 equation^{22–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 Lauritsen^{27} 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 $\Phi $, field enhancement $\beta $, 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 Eyges^{31}) 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 ($\beta $) and apex radius ($an$) factors. Using the PCM to construct a hyperbolic or ellipsoidal parameterization enables utilizing findings from prolate spheroidal models^{11,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 photocathodes^{15,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 model^{21} is improved and applied to treating randomized surface roughness and allows roughness contributions to emittance for photocathodes to be evaluated.

## II. ARRAY MODEL

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 $12CaO\u22127Al2O3$:$4e\u2212$, 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\u2212$, 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 $\xb0C$ 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 $\u223c$$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.

### A. Emitter and array geometry

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\mu m$, a base diameter of $2.60\mu 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\xd7$ 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\xd7$ via Schottky’s conjecture^{13,36,37}) is not included for nanocages on the conical surface.

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\mu m$ pitch (emitter density $1.28\xd7107cm\u22122$) would require 78 nA per identical emitter to produce a $1A/cm2$ current density. For typical temperatures of interest in heated cathodes ($650\xb0C\u2013900\xb0C$), 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.

In the PCM approach, a single emitter is composed of $n$ point charges in a dipole arrangement^{37,54,55} where the charges are placed along the $z^$-axis at positions $zj$ ($j\u2208[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 $Fo\u2261q|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 $x$–$y$ plane^{56} 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 $\rho ~$ and $z~$ are the dimensioned coordinates, then let $\rho ~=\rho lo$ and $z~=zlo$, where $\rho $ and $z$ are understood to be dimensionless. Lastly, the locations of the point charges are placed such that

which implies that $zj=zj\u22121+rj$ and that the height of an isolated emitter is

The factor $r$ governs the sharpness and height of the emitter by dictating the relative magnitudes of the $n$ point charges^{57} (see Appendix B). The energy scale $Folo$ is introduced. The dimensionless function $Vn(\rho ,z)$ is then defined by

where $U$ is the dimensioned potential energy, the radius factors $rj\xb1$ are distinct from the ratio factor $r$, and $\lambda j(r)$ is proportional to the magnitude of the charge at $zj$. $\lambda $ absorbs factors such as $q/4\pi \epsilon 0$ and others that otherwise would appear: it vanishes rapidly as $j$ increases for small $r$, and, therefore, defining $\lambda j(r)\u2261rjPj(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(\rho ,z)$ is analytic, and the apex radius $an$ and field enhancement $\beta 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 $\beta $ is the ratio of the local field $F$ on the surface [defined by $zs(\rho )$] to the background field $Fo$, and it is, therefore, dimensionless and may be calculated via gradients of $Vn$ with $z$ and $\rho $. Observe that $an$ is dimensionless so that $anlo$ is the dimensioned apex radius. It exponentially decreases with $n$ for $r<1$ just as $\beta n$ exponentially increases. Matrix methods for the evaluation of $\beta n(r)$ and $an(r)$ for the apex of an isolated dipole PCM emitter are summarized in Appendix B. As $r\u21921$, the shape of the resulting emitters approaches the shape of an ellipsoid: although $r>1$ solutions are possible, the line charge model^{58} 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 base^{54} 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\Delta ,y0+j\Delta )$. 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 $\u27e8Itip\u27e9$ and the array density of emitters $\rho array$ such that $Jarray=\u27e8Itip\u27e9\rho array$. For the same array density, a triangular ($t$) and square ($s$) array will have different pitches given by

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 $\Delta \xd7\Delta $. In contrast, a triangular lattice is equivalent to two overlapping rectangular grids with grid elements of size $\Delta \xd73\Delta $ with the origin of the first rectangular lattice at $(0,0)$ and of the second at $(\Delta /2,3\Delta /2$).

Consider a circular macroscopic emission area with the $x$-axis along the horizontal diameter. Let $nt\u22612k+1$ be the number of emitters along that diameter, spaced every $\Delta $ (that is, $xi=i\Delta $ for $\u2212k\u2264i\u2264k$). For a square ($s$) lattice, the placement of lattice sites above and below the diameter is governed by $yj=j\Delta $ for $\u2212nt2\u2212i2\u2264j\u2264nt2\u2212i2$ such that the diameter of the total emission area is $nt\Delta =(2k+1)\Delta $. 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^(\varphi )$ and a basis vector $b\u2192ij$ defined by

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

An initial rotation angle $\varphi ij$ is defined by

A hexagon of sites is then generated by initially rotating $b\u2192ij$ by $M^(\varphi ij)$ followed by the application of $M^(\pi /3)$ six times to find the sites of the six points of a hexagon. That is, the site of the $l$th vertex of the hexagon is located at $M^(\pi /3)l\u22c5M^(\varphi ij)\u22c5b\u2192ij$. 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 $k\u22646$, these are the same as arrays allowing all sites within a circular region of diameter $(2k+1)$. For $k\u22657$, 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 $\beta $ and apex radius $a$ of an emitter in our array model changes depending on pitch $\Delta $ (an isolated emitter is equivalent to $\Delta \u2192\u221e$). 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, $\beta 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, $\beta (\rho )$ will refer to the local field enhancement along the surface of the rotationally symmetric emitter defined by $zs(\rho )$, where $\rho 2=x2+y2$, and $a$ will denote the apex radius of the curvature for general pitch $\Delta $. At the apex of the $jth$ emitter, they are symbolically given by (generalizing Ref. 55)

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)\u2192Vn(\rho ,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 approximations^{60} so that if $fj=U(x,y,z+j\delta )$ with $\delta \u226a1$, then the $kth$ iteration to find the surface is given by

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, $\beta (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.

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 models^{11,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.

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 $\rho \u2192x$), let $z1=zs(0)$ and $z2=zs(x2)$. The pairs $[xj,zj=zs(xj)]$ define the numerically determined points of the surface and $x1\u22610$ 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

in terms of which

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

The field enhancement factor $\beta (x)$ as a function of position $x$ along the surface for a hyperbolic emitter is analytic ($\rho \u2192x$ for $y\u21920$). That it is a reasonable approximation to the PCM $\beta (x)$ is demonstrated as follows. The hyperbolic emitter field enhancement factor^{62} is given by

where $\upsilon 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\eta sin\upsilon o$, then

where $\beta j=\beta (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.65\u2272r\u22720.85$ for $n\u22485$. 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.,

with $sinh\eta 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 $\Delta $ 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 $\Delta =5$ $\mu $m entails $k\u22481000$. Emitters around the rim have negligible area of contribution: the outer edge of emitters for a $D=1$ cm diameter disk accounts for $(\Delta /D)(2\u2212\Delta /D)<0.01$% of the emission area. A large area array can, therefore, be taken as the $k\u2192\u221e$ limit with respect to shielding effects. The $\beta k\u2192\u221e$ limit may be estimated by observing that $\beta k/\beta 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 $\Delta =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=\Delta /Sn+1(r)$ for $n=5$ and $r=0.7$ [for which $S6(0.7)=2.9412$], shows the diminishing effect of shielding as $\Delta $ 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 ($p\u22656$) is within 4% of the unshielded case (close to the $p=12$ value).

$p$ . | Sq Car . | Hex Car . | Sq Tri . | Hex Tri . |
---|---|---|---|---|

2 | 0.309 553 | 0.305 726 | 0.259 846 | 0.259 088 |

3 | 0.559 134 | 0.557 655 | 0.526 304 | 0.525 689 |

4 | 0.657 410 | 0.656 726 | 0.639 896 | 0.639 570 |

5 | 0.698 481 | 0.698 116 | 0.688 609 | 0.688 429 |

6 | 0.717 969 | 0.717 752 | 0.711 976 | 0.711 869 |

7 | 0.728 255 | 0.728 117 | 0.724 379 | 0.724 310 |

8 | 0.734 156 | 0.734 063 | 0.731 516 | 0.731 470 |

9 | 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 |

$p$ . | Sq Car . | Hex Car . | Sq Tri . | Hex Tri . |
---|---|---|---|---|

2 | 0.309 553 | 0.305 726 | 0.259 846 | 0.259 088 |

3 | 0.559 134 | 0.557 655 | 0.526 304 | 0.525 689 |

4 | 0.657 410 | 0.656 726 | 0.639 896 | 0.639 570 |

5 | 0.698 481 | 0.698 116 | 0.688 609 | 0.688 429 |

6 | 0.717 969 | 0.717 752 | 0.711 976 | 0.711 869 |

7 | 0.728 255 | 0.728 117 | 0.724 379 | 0.724 310 |

8 | 0.734 156 | 0.734 063 | 0.731 516 | 0.731 470 |

9 | 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 $\Delta $ (in contrast to prior studies using the LCM^{1,29,58} where the apex radius was held fixed): as the equipotential curves become progressively less blunt as $\Delta $ decreases, $\beta k$ (which is proportional to the inverse apex radius) correspondingly decreases. As an example, for $\Delta =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 $\beta k$ with decreasing $\Delta $ would be less.

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 $\lambda 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.

### B. Emission parameters

The behavior of $\beta (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:

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;

when the pitch $\Delta $ 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 ($\Phi \u2192\Phi \u22124QFo$ in the notation of Ref. 24) in the RLD equation overestimates the current;

inferring $\Phi $ 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,

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

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 $\delta $ is situated at the center of a circle of diameter $\Delta $. The current density is

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, $JGTF\u2192JFN$.

A region adjacent to the apex is characterized by $n=1$ for which $N(1,s)=(s+1)e\u2212s$.

Current density near the base of the emitters is characterized by $n<1$ for which $JGTF\u2192JRLD$.

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)\Delta 2\u2212\pi \delta 2$.

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

At $n=1$, $No(1,s)=s$ exactly because the vanishing of $e\u2212(n\u22121)s\u22121$ offsets the singularity in $\Sigma (n)$ as $n\u21921$. It is seen that

which exceeds unity. The accuracy of $\Delta N$, therefore, determines the accuracy of Eq. (A19). However, $s$ is large: when $n=1$, then $s\u224840$ for metal-like work functions ($\Phi =4.5eV$) and $s\u224810$ for low work functions ($\Phi =2eV$) when $F=1eV/nm$. As a result, the small error $(\Delta N\u22121)/s\u226a1$ 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 $n\u21921$ 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 $\rho j$ for the $jth$ ribbon. There is no requirement that $\rho j$ be evenly spaced: near the apex of the emitter where fields vary most strongly, $\rho j$ is more densely spaced than further from the axis of symmetry. Field enhancement values $\beta j=\beta (\rho j)$ may be constructed from Eq. (15) with the replacement $x\u2192\rho $. The current density at the locations $\rho j$ is then defined by the components

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

For a surface $zs(\rho )$ specified by the discrete points $(\rho j,zj)$ with $zj=zs(\rho j)$, and approximating the gradient of $zs(\rho )$ using finite differences,^{60} then

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

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

The usage of the hyperbolic surface $zh(\rho )$ simplifies the model but introduces an ambiguity indicated in Fig. 9: at some point, $\beta h<\beta (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 $\beta (\rho )\u2192\beta (\Delta /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

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 $\beta (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.

Parameter . | Symbol . | Value . | Unit . |
---|---|---|---|

Temperature | $T$ | 1173 | K |

Work function | $\Phi $ | 2.3 | eV |

Fermi energy | $\mu $ | 7 | eV |

Length unit | $lo$ | 1 | $\mu m$ |

Charge number | $n$ | 5 | $\u2026$ |

Self-sim ratio | $r$ | 0.7 | $\u2026$ |

Pitch | $\Delta $ | 8.8235 | $\u2026$ |

Tip number | $k$ | 3 | $\u2026$ |

Height | $zs(0)$ | 2.9412 | $\u2026$ |

Base radius | $\delta $ | 1.7071 | $\u2026$ |

Lattice | $\u2026$ | Tri | $\u2026$ |

Area | $\u2026$ | Hex | $\u2026$ |

Parameter . | Symbol . | Value . | Unit . |
---|---|---|---|

Temperature | $T$ | 1173 | K |

Work function | $\Phi $ | 2.3 | eV |

Fermi energy | $\mu $ | 7 | eV |

Length unit | $lo$ | 1 | $\mu m$ |

Charge number | $n$ | 5 | $\u2026$ |

Self-sim ratio | $r$ | 0.7 | $\u2026$ |

Pitch | $\Delta $ | 8.8235 | $\u2026$ |

Tip number | $k$ | 3 | $\u2026$ |

Height | $zs(0)$ | 2.9412 | $\u2026$ |

Base radius | $\delta $ | 1.7071 | $\u2026$ |

Lattice | $\u2026$ | Tri | $\u2026$ |

Area | $\u2026$ | Hex | $\u2026$ |

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(\Delta /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 $\phi a$ (e.g., $Fo=q\phi a/D$ for which $\phi a=2kV$ if $D=100\mu 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 area^{62} of $\pi as2g[\beta (0)Fo]$ where for hyperbolic emitters,^{23,64}

where^{65} $b(F,\Phi )=42m\Phi 3/(3\u210fF)$ and $\nu (\Phi )=(8Q/9\u210f)2m/\Phi $. For a cone angle of $tan\u2061(\theta )=1.3/1.6$ (or $\theta =39\xb0$) and $\Phi =2.3eV$, $g(F)$ simplifies to $g(F)=0.88150F/(34.869\u2212F)$ for $F$ in units of eV/nm. Using the RLD equation [Eq. (A20)] instead of $JGTF$, then for the same total area characterized by $\Phi $ and a surface field of $\beta (\Delta /2)Fo$ where $\beta (\Delta /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=\beta (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\u2061(I)$ vs $F$], and in FN coordinates [$ln\u2061(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.

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 $\Phi eff$ is determined from the slope $\u2212\Phi eff/kB$ on an RLD plot of $ln\u2061(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 $\Phi =2.3eV$, the gray lines would result and the slopes would correspond to $\varphi =\Phi \u22124QFo$. With the field enhancement at the apex of $\beta =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.

For field data, the effective field enhancement factor $\beta eff$ is determined (if $\Phi $ is known) from the slope on an FN plot of $ln\u2061(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 $\u2212b(Fo,\Phi )/\beta $, where $b(F,\Phi )$ is as encountered in Eq. (25) and the variation of the other factors in $IFN$ with $\beta $ 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(\beta Fo)$, the single gray line would result and the slope would correspond to $\beta =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.

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.

## III. ROUGHNESS MODEL

### A. Random surface specification

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 $\beta $-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

so that the simulation region is centered at $(x,y)=(0,0)$ and spans $\u22121/2\u2264(x,y)\u22641/2$, and the outermost point charges are a distance $\delta $ 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 $\u2212zo$ (dipole model). To randomize the locations, a radius and an angle characterized by $\eta 1$ and $2\pi \eta 2$ are added such that $xi\u2192xi+\eta 1cos(2\pi \eta 2)$ and $yj\u2192yj+\eta 1sin(2\pi \eta 2)$, where the $\eta j$ are uniformly distributed random numbers $0\u2264\eta j\u22641$. The potential energy is then

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

where

where $U$ is the (dimensioned) potential energy, and the summations over $k$ and $l$ range from $\u2212n$ 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 $\u2202zV\u2248[V(z+\u03f5)\u2212V(z\u2212\u03f5)]/2\u03f5$ so that the equation $V(z+\u03f5)=0=V(z)+\u03f5\u2202zV$ is iteratively solved by

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.

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

where $xj$ is $x$, $y$, or $zs(x,y)$. The total force $F2=Fx2+Fy2+Fz2$ is the product of a field enhancement factor $\beta (x,y,zs)$ with the background field $Fo$ such that $|F\u2192|=\beta Fo$. The variation of $\beta (x,y,zs)$ for the surface of Fig. 16 is shown in Fig. 17, and it is used to find the local field $|F\u2192|=\beta Fo$ used in the emission equation $JGTFP(F,T)$. The enhancement is seen to be above background ($\beta >1$) near the apexes of the protrusions, but significantly, “suppression” ($\beta <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}

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

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)\delta 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

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

on each element with $\Phi =2.3$ eV, $T=300$ K and $Fo=0.05$ eV/nm, $4QFo=0.268325$ eV, and $\lambda =532$ nm ($\u210f\omega =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.

### B. Trajectories

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\rho $ 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\u2212(za\u2212zc)(j/N)$ as intermediate locations: the variation of the *rms* departure as a function of $j$ is then given by

where $Vj\u2261V(x,y,zj)$ and the mean $\u27e8Vj\u27e9$ 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 $\u27e8Vj\u27e9$ by $z\u22482zc$. 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.

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

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

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

Evaluate $\beta (x,y)$ at $zs$.

Specify $\sigma $ and $\kappa $ via Eqs. (37) and (D8), and $C$ by comparisons to trajectories [or Eq. (D18)].

Evaluate $v\u2192(t)$ for launches distributed over surface.

Specify the time $ta$, defined by $z(ta)=za\u223c2zc$.

Weight $v\u2192(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)\u2248\sigma (x,y)v\u2192o\u22c5k^$ for the $k^th$ component, which allows $vz(ta)$ in the Ballistic-Impulse approximation to be given by Eq. (D17), rewritten as

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^\u22c5z^\u2261cos\theta =Fz/\beta Fo$, and so, $vz(0)=voFz/\beta Fo=vocos\theta $ where^{71} $FzL=\u2202zU|z=zs$, then $vx(0)=vosin\theta cos\varphi $ and similarly for $vy$, to give

where $\varphi $, like $\theta $, 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

if the current density does not change appreciably over the time $\Delta 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 $\beta F\u226b1$ 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

where $\u230a\cdots \u230b$ 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 $\u27e8vk2\u27e9$ (from which “mean transverse energy” factors are evaluated) is then

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 $v\u2192\u22c5k^$ defined in Eq. (43) are obtained as follows. Define $zs(xi,yj)\u2261zij$ on a grid such that $xi=(i\u22121)\Delta $, $yj=(j\u22121)\Delta $, and $\Delta =L/(Np\u22121)$ and $Np$ the number of grid points along an axis. Introduce two diagonal vectors,

The surface normal vector $n^$ for $Aij$ is then

[not to be confused with the weighting factor of $nij$ of Eq. (42)] from which the angles $\theta ij$ and $\varphi ij$ are defined by $cos\theta ij\u2261n^ij\u22c5z^$ and $cos\varphi ij=(n^ij\u22c5x^)/sin\theta $, the magnitudes of which are

Alternately, if $(1/2)mv(ta)2$ is the kinetic energy at the anode plane, then the $x^$-component of the $v\u2192a=v\u2192(ta)$ is $v\u2192a\u22c5n^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 $\kappa =72.40$ (parameters characteristic of photoemission as given in Table III of Ref. 21).

Histograms of the transverse velocity component $v\rho 2\u2261vx2+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 $\alpha $ such that the normalized distribution $f(v\rho )$ resembles

where $\alpha $ may be crudely estimated by the location $v\rho =vmax$ where $f(vmax)$ is the maximum value of the numerically evaluated ballistic-impulse histogram or $\alpha =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\rho )$ may be used to estimate the MTE, using the trajectory weightings and velocity components directly to obtain

is preferable.

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 $MTE\u223c\u27e8Et\u27e9$ 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 Schmerge^{72} that experimentally measured thermal emittance is approximately a factor of two greater than the theoretical estimate.

### C. Comparison to particle-in-cell

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 MICHELLE^{73} is used to validate the SRM theory. The current density weighting factor $w$, field enhancement factor $\beta $, and transverse velocity distribution $v\rho $ 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.

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.

The field enhancement factor $\beta $ 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.

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.

## IV. CONCLUSION

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.

## ACKNOWLEDGMENTS

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.

### APPENDIX A: GTFP CURRENT DENSITY

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 treatments^{22–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 ($\beta T$) and field ($\beta F$) slope factors are obtained from the gradients with energy of the exponents of the Fermi-Dirac distribution (thermal) and the Gamow factor $\theta (E)$ on which transmission probability is based (field). The energy slope factor $\beta T$ for temperature is independent of energy and given by

where $\mu $ is the Fermi energy. The energy slope factor $\beta F(E)$ for field depends on energy and is related to the gradient of the Gamow factor $\theta (E)$ by

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

where $U(x)$ is the one-dimensional potential energy including the image charge modification, $U(x\xb1)\u2212E=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

where $F$ and $Q$ are defined by these equations. The Gamow factor $\theta (E)$ is taken as a cubic in energy for energies above the Fermi energy and below the barrier maximum,^{22} that is, for when $\mu \u2264E\u2264\mu +\Phi \u22124QF$, which implies that $\beta F(E)$ is a quadratic in energy. The coefficients of the interpolating polynomial can be inferred from the behavior of $\theta (E)$ near $E=\mu $ (which is well approximated by the Fowler-Nordheim representation) and near $E=\mu +\Phi \u22124QF$ (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 integrals^{26} but for which very good approximations are^{75}

where $y\u22614QF/\Phi $ and is related to the Schottky-reduced barrier term,

Coefficients of the cubic polynomial can be uniquely determined by specifying the values $\theta (\mu )$ and $\theta (\mu +\varphi )$ and derivatives $\beta F(\mu )$ and $\beta F(\mu +\varphi )$ at the end points, which are known to be^{62,76}

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

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

For $T<Tmin$:

set $n=\beta T/\beta F(\mu )$ and $s=\theta (\mu )$

For $T>Tmax$:

set $n=\beta T/\beta F(\mu +\varphi )$ and $s=\beta F(\mu +\varphi )\varphi $

Otherwise:

set $n=1$, let $Em$ be defined by $\beta T=\beta F(Em)$ and(A13)$s=\beta T(Em+\theta (Em)\beta F(Em)\u2212\mu ).$

Construct

where $ARLD=(qmkB2)/(2\pi 2\u210f3)=120.1735A/cm2K2$. Let $N(n,s)$ be defined by

It can be shown that

where $\zeta (2)=\pi 2/6$ (Riemann Zeta function). A useful approximation is enabled by introducing

in terms of which

The following limits represent the connection with the canonical equations:^{23}

which is tantamount to $n\u21920$ and $n\u2192\u221e$, respectively. Equation (A20) is the Richardson-Laue-Dushman (RLD) equation^{25} when small reflection coefficients are neglected and an “effective” work function $\varphi =\Phi \u22124QF$ is used,^{77} and Eq. (A21) is the Murphy and Good^{26} 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 $\Phi $ and $\beta $ values derived from $I(V)$ data.

$F$ . | $T$ . | $\Phi $ . | $n$ . | $s$ . | $log\u2061JGTF$ . | $log\u2061JRLD$ . | $log\u2061JFN$ . |
---|---|---|---|---|---|---|---|

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 |

$F$ . | $T$ . | $\Phi $ . | $n$ . | $s$ . | $log\u2061JGTF$ . | $log\u2061JRLD$ . | $log\u2061JFN$ . |
---|---|---|---|---|---|---|---|

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 $(\mu +\varphi \u2212E)$ appearing in the transmission probability has $E$ augmented by the photon energy $\u210f\omega $ which renders $(\mu +\varphi \u2212E\u2212\u210f\omega )<0$. Under typical photoemission conditions, $N(n,s)$ now becomes negligible compared to $N(n,\u2212s)$, and so the emitted current density is proportional to^{23,24}

### APPENDIX B: FIELD ENHANCEMENT AND RADIUS

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\u2032$ and $P$ for the dipole defined by^{55}

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

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 $\beta n$ is found from

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

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, $\beta n/\beta N\u22481.24n\u2212N$, reflecting Schottky’s conjecture regarding the multiplicity of $\beta $-factors.^{37} Similarly, $aN/an\u22480.72N\u2212n$ for the chosen value of $r$. The exponential behavior of both is shown in Fig. 30.

$n$ . | $Sn(r)$ . | $Pn(r)$ . | $\beta n(r)$ . | $\beta n/\beta n\u22121$ . | $an(r)$ . |
---|---|---|---|---|---|

1 | 1.0000 | 2.295 0 | 4.0582 | $\u2026$ | 0.664 53 |

2 | 1.7000 | 1.537 2 | 5.0640 | 1.2479 | 0.557 21 |

3 | 2.1900 | 1.308 5 | 6.3931 | 1.2625 | 0.420 66 |

4 | 2.5330 | 1.131 2 | 8.0182 | 1.2542 | 0.308 62 |

5 | 2.7731 | 0.979 60 | 9.9999 | 1.2472 | 0.223 29 |

6 | 2.9412 | 0.847 85 | 12.424 | 1.2424 | 0.160 19 |

7 | 3.0588 | 0.733 29 | 15.398 | 1.2393 | 0.114 27 |

8 | 3.1412 | 0.633 87 | 19.053 | 1.2374 | 0.081 171 |

9 | 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 |

$n$ . | $Sn(r)$ . | $Pn(r)$ . | $\beta n(r)$ . | $\beta n/\beta n\u22121$ . | $an(r)$ . |
---|---|---|---|---|---|

1 | 1.0000 | 2.295 0 | 4.0582 | $\u2026$ | 0.664 53 |

2 | 1.7000 | 1.537 2 | 5.0640 | 1.2479 | 0.557 21 |

3 | 2.1900 | 1.308 5 | 6.3931 | 1.2625 | 0.420 66 |

4 | 2.5330 | 1.131 2 | 8.0182 | 1.2542 | 0.308 62 |

5 | 2.7731 | 0.979 60 | 9.9999 | 1.2472 | 0.223 29 |

6 | 2.9412 | 0.847 85 | 12.424 | 1.2424 | 0.160 19 |

7 | 3.0588 | 0.733 29 | 15.398 | 1.2393 | 0.114 27 |

8 | 3.1412 | 0.633 87 | 19.053 | 1.2374 | 0.081 171 |

9 | 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 |

### APPENDIX C: SPHERICAL AND PROLATE SPHEROIDAL MODELS

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

For the prolate spheroidal coordinate system,^{62} then

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

The spherical case is recovered in the limit $\upsilon \u2192\u221e$ for which $ae/be\u21921$.

The hyperbolic coordinates are defined by $(z/bh)2\u2212(\rho /ah)2=1$ with $ah=Lsin\upsilon $ and $bh=Lcos\upsilon $ in Eq. (C2) so that

### APPENDIX D: BALLISTIC-IMPULSE MODEL

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,\theta )$ with boundary conditions such $U(a,\theta )=0$, giving

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

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,\theta )=\beta (\theta )Fo$ is then seen to be $\beta (\theta )=3cos\theta $ 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=\rho 2+z2$ and $tan\u2061\theta =\rho /z$ (where $z$ and $\rho $ have dimensions of length while finding the ballistic equations) are given by

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 $\rho =0$ and $z/a=(2,3)$, $Fz/Fo=(1.25,1.07)$.

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

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

The ballistic equations are then tantamount to $F\rho \u21920$ and $Fz\u2192Fo$, with the launch sites defined by $z(0)=acos\theta $ and $\rho (0)=asin\theta $. Introduce $p\u2261vot/a$ as a dimensionless time and $\kappa $ as the ratio of the energy scale associated with the background field and that associated with the initial kinetic energy defined by

where $\tau =a/vo$, then the ballistic trajectories are defined by $(\rho b(t),zb(t))$, and are given by

The “Impulse Approximation” augments the initial velocity $vo\u2192\sigma 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

Assuming the impulse occurs only close to the hemisphere for a time $p=a/vot\u223c1$ identifies the impulse factor $\sigma $ as

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 ($\theta =0$), a particle passing the location $z=na$ for $n\u22482$ 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

Substituting the first into the second equation, letting $\sigma \u22481+(3/2)C\kappa $ (on axis), assuming $n2\u226b1$, and collecting gives

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.

## REFERENCES

*Materials Research Society Symposium Proceedings*(Materials Research Society, Warrendale, PA, 1998), Vol. 498, pp. 197–208.

*IEEE International Vacuum Electronics Conference (IVEC)*(IEEE, 2018), pp. 39–40.

*Proceedings of FEL Conference*(BESSY, Berlin, Germany, 2006), pp. 583–586.

*International Particle Accelerator Conference*(IEEE, New Orleans, LA, 2012), p. MOPPP041.

*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).

*Wiley Encyclopedia of Electrical and Electronics Engineering*, edited by J. G. Webster (John Wiley & Sons, Inc., New York, 2014).

**2**, 123–242 (1930) who give the RLD equation [their Eq. (6)] in the form we have, identify $\varphi =\Phi \u22124QF$ as the “effective” work function as we have, and use $F=q|E|$ as we do.