Thermal-field emission from cones and wires

The generation and manipulation of high brightness electron beams is directly dependent upon the physics of the emission process, but beam simulation presumes instantaneous emission using 1D semianalytic theories applied to simplified static barriers. The prediction of emission and beam properties is essential for the predictive design and development of technologies of critical importance to devices that rely on tunneling and/or vacuum transport, particularly in nanodevices^1,2 and harsh environments, such as pulses drawn from carbon fibers^3,4 and hollow cathodes for charge neutralization in satellites.⁵ More generally, emission modeling is used to simulate many technologies relevant to commercial, military, and research needs: High Power Microwave devices,⁴ directed energy,³ RF/THz technologies,^6,7 electronic warfare, multi-$e$beam lithography and electron microscopy,^8–11 radiation-hard and high temperature electronics, energy harvesting,^12,13 rf and voltage breakdown,^14,15 microdischarges and the generation or suppression of cathodic arcs,¹⁶ ion propulsion/microthrusters for satellites,^5,17,18 and satellite charge neutralization¹⁹ and electrodynamic tether applications.^20,21

Field emission enables high brightness electron sources that are either applied to or under development for applications central to research, industrial applications, and military technologies. In some applications, a combination of temperature and curvature effects makes the characterization and modeling of these sources both important and comparatively difficult. Generally, the electron emission equations are taken to be one of the canonical mechanisms, but the ignoring of temperature and curvature effects as generally done undercuts simulation accuracy and hence corrupts presumed beam characteristics. Such concerns are not problematic for relatively blunt (or flat) emitters run either very hot or at current levels where heat generation is not an issue. Nevertheless, for high brightness sources, in particular, where high current density and narrow angular spread of the beam are desired, multiple effects conspire to make their combined modeling difficult even if their underlying complications are individually well accounted for.

When the emitter is or becomes hot, thermal-field emission sources that can meet taxing performance demands, withstand radiation-hard and poor vacuum environments, and enable thermally tolerant electronics will contribute to many emerging technologies.²² Performance, ruggedness, and survivability are stressed. For nanogap devices [e.g., nanovacuum channel transistors (NVCTs) and photo-excited structures^1,2,23,24], a vacuum channel circumvents technological obstacles caused by, first, significant Joule heating in the transport channel and Nottingham heating at the emission site²⁵ and, second, radiation damage to electronic materials. Ballistic transport in a vacuum channel does not generate heat, and realizable vacuums have an enormous breakdown threshold (⁠$Ebd∼100$ MV/m). The ratio of the anode-cathode (AK) separation $W$ with the transit time $τ$ replaces the semiconductor saturation velocity $vs→W/τ$ in Johnson’s figure-of-merit²⁶ (⁠$J=Ebdvs/2π$⁠): for example, NVCTs enable $Jvac>103×JSi$ (in contrast, $JGaN=15.4×JSi$⁠). Emission into vacuum is by quantum mechanical tunneling through a surface barrier (⁠$Φ$⁠), but the physics applies equally well to microscopy, heterostructure, nanodiode, and metal-insulator-metal (MIM) structures.^22,27–29 Tunneling occurs when the DeBroglie wavelength $λB$ of an electron (⁠$−q$⁠) satisfies $qE≈Φ/λB∼9$ eV/nm. Nanodimensioned AK gaps (smaller than mean free path) extend advantages by substantially relaxing vacuum requirements. For nanofiber and nanotubes, heating changes field to thermal-field emission^30,31 because of strong temperature gradients and nonlinear heat loss factors, making prediction very difficult and undermining characterization of emission properties using conventional steady state models.^25,32 Dynamic effects of ions approaching surfaces influence breakdown behavior in poor vacuum (Paschen’s law). Clearly, in many circumstances, not all thermal-field emission is an opportunity: undesirable microprotrusions on metal surfaces contribute to high-field breakdown in high power devices and is described by similar physics.

The present study describes efforts to create theoretical models of thermal-field emission with position-dependent temperature variation and dynamic processes using methods derived from quantum distribution, transfer matrix (tunneling), thermal-field emission, and heat transport/radiative loss methods, so as to develop computationally rapid algorithms appropriate for beam optics codes treating ultrafast/ultrasmall processes. The models are to address the simulation of ultrashort electron pulses in vacuum for nanodevices and rapidly heated nanotube/fiber emitters that have complex emission barriers and dynamic launch conditions associated with time-dependent and nonuniformly heated surfaces. Nanodevices entail novel tunneling barriers affected by charge in the AK gap, which leads to “noise” and stochastic processes. Nanofibers have temperature gradients along the fiber that can approach “runaway” conditions.²⁵ The apex of fibers can taper,³³ analogous to a hyperbolic emitter,³⁴ and can, therefore, require a locally varying image charge contribution to the barrier.^35–37 Many time scales are involved: tunneling time, ultrafast field variation, electron relaxation time during transport, heat diffusion times from Joule and Nottingham effects but opposed by radiative heat loss processes, and nanoscale AK gap transit times. They depend on current governed by local thermal-field (TF) emission processes^30,31 past barriers associated with coatings and heterostructures²² and space-charge oscillations. Computationally efficacious models are desired for, e.g., particle-in-cell^38–40 and molecular dynamics^41,42 codes used to describe electron sources with strong TF and temporal variations as part of their normal operation in a manner allowing for treating the consequences of finite tunneling, scattering, and AK gap transit times on the emission process and result in intensive efforts to speed up computations associated with local and variable emission conditions.^32,37

Emission events in simulation codes are commonly treated as instantaneous (applied fields are static) because the relevant time scales of changes are generally long compared to those of emission itself (tunneling, scattering, transit times, radiative processes): simulations of wave packets upon barriers using a quantum distribution function (specifically, the Wigner distribution function) show that the packets are delayed in time compared to trajectories that follow a ballistic model. More egregiously, current at the presumptive emission site is treated using only one of the asymptotic current limits associated with either field or thermal emission, when, in fact, emission from thermal-field conditions can be orders of magnitude larger, and past barriers that depart from the assumed “image charge potential barrier” of the canonical emission equations. These complications attend nanogaps^43,44 and short bunches (particularly when from nanotubes and fibers),³² in particular, for which the static environment/instantaneous emission/field or thermal approximations are ill suited. Recent developments using the Schrödinger equation^22,31,45 and quantum distribution function^46,47 methods provide a basis to formulate thermal-field and time-dependent emission models, which allow for the simulation of dynamic processes, e.g., ones that can contribute to breakdown and/or thermal runaway along fibers,²⁵ and offer reasonable models of delays due to overbarrier and tunneling transport.

Relevant time scales are ultrafast bunches (10–1000 fs), scattering in bulk (1–100 fs), electron AK gap transit times⁴³ (electrons travel 100 nm in 15.1 fs for $E=5$ GV/m), and Büttiker–Landauer semiclassical times (⁠$τsc∼Φ/qEvF≈0.6$ fs) characterize tunneling delay. Noise (random emission events), intrinsic emittance (reduces beam brightness), and launch angle (causes gate interception or beam spread) depend on temperature and the local variation of the emission barrier due to coatings and heterostructures.¹⁸ To connect dynamic quantum behavior to trajectory and beam optics, emission models must provide

• the dependence of current density on the local field and temperature conditions, which change along the structure,

• the resulting field emission velocity distribution from the surface,

• the mean transverse energy at the emission site, which is the dominant component of intrinsic emittance, and

• the launch time of the particle once chosen to be emitted.

All manifest as $I(V)$ hysteresis, duty cycle variations, or time-dependent evolution, and therefore, conventional asymptotic steady state derived algorithms in simulation are undermined.¹⁸ 

Moreover, semiconductor field emission (as from carbon fibers) entails significant heating that can occur along fibers, nanotubes, or high aspect ratio nanowires when driven to large per-tip currents^25,31 such that existing models cannot anticipate thermal runaway (emitters are commonly current-limited to a 1 $μ$A/tip to prevent failure in nanogap devices). Variations in doping and quantum mechanical confinement (discrete energy levels associated with triangular wells depend on doping) change band bending effects at the surface, and the presence of oxides alter the tunneling barrier and make the tunneling theory substantially more difficult.^22,45 Dynamic theoretical emission models must be computationally tractable to be useful to device and beam simulation: density functional theory is far too numerically intensive, and Transfer Matrix Approach (TMA) methods, even though fast, are still too time-consuming. In simulation, a surface mesh describing just the emitter has millions of elements contributing to total current in a simulation, each of which contains different conditions for the emission model. The proposed theory herein enables predictive simulations of ultrafast and ultrasmall emitter conditions and trajectory mapping by creating parametric and computationally rapid theoretical models vetted by more comprehensive theoretical methods. Relating emission to dynamic processes changes the modeling paradigm significantly. A program with several long term aims is recommended to develop emission velocity and launch time models. Doing so would meet several needs:

• specification of launch conditions for simulations that account for dynamic changes in the emission barrier due to field variation and emission fluctuation,

• evaluation of self-consistent temperature variation along nanowires caused by the coupled general thermal-field current (GTF),³⁰ 

• consideration of nonstandard barrier models vetted using TMA methods, and

• accounting for launch delay associated with passage over and through barriers.

Neither the TMA method nor even exact evaluations of $θo(E)$ or its generalizations below are suitable for inclusion in modern high performance particle-in-cell (PIC) codes^40,48–51 or molecular dynamics^41,42 because the computational cost of numerical evaluations makes them unattractive, if not prohibitive, for simulations requiring emission from millions of surface elements over millions of time steps in time-dependent and spatially nonuniform simulations. Thermal-field emission frustrates the use of purely thermal or field emission models because $θo(E)$ departs from linear approximations near the Fermi level or barrier maximum. The GTF equation and its refinements^30,52 were designed to address both computational expediency and accuracy particularly when thermal-field processes contribute^18,31 (see Ref. 37 for a complementary approach that provides greater speed). In these cases, however, the Schottky–Nordheim (SN) barrier is presumed. The need remains to develop a high speed, high accuracy method analogous to the GTF method for non-SN barriers (particularly so-called quadratic barriers) and to provide an assessment of predicted consequences for when that model is applied to a more realistic model of a fiber or nanotube field emitter where the apex is not cylindrically flat but rather tapered to a point, and curvature plus space-charge forces alter the barrier. The present study will address those needs by

• using accurate TMA methods to construct generalized $θ(E)$ factors,

• developing an accurate high speed method based on the shape factor approach to evaluate conventional $θo(E)$ factors,

• extending the shape factor method to generalized $θ(E)$⁠, particularly for quadratic barriers,

• assessing the impact on TF emission by considering the degree to which temperature increases for a conical emitter beyond the cylindrical model, and

• indicating the probable consequences for the total thermal-field current from such an emitter.

The TMA method for evaluating the transmission probabilities uses the methods for matching wave functions constructed from Airy function solutions and first derivatives at locations where a discontinuous change in barrier height, a change in the slope of adjacent barrier segments, or both is encountered.^33,53,54 It accurately and exactly accounts for both tunneling and fly-over contributions to the emitted current and, therefore, accounts for both field and thermal emission contributions. It is applied to representative tunneling barriers and then focused, in particular, on the standard Schottky–Nordheim barrier (also known as the image charge barrier) for metal-like parameters to address the adequacy of the Gamow factor and its use in the Kemble approximation for the transmission probability and their usage in the GTF equation.³⁰ For simplicity, energy will always be taken as parabolic in momentum or $E(k)=ℏ2k2/2m$ in the one-dimensional tunneling equations. Units throughout follow the conventions of prior studies, particularly, Refs. 33 and 54, and, in particular, use units such that energy is in eV, distance is in nm, time is in fs, and the unit (positive) charge $q$ is absorbed into the definitions $F≡q|E|$ and $V(x)=qφ(x)$⁠, where $(E,φ)$ are the electric field and potential, respectively. As a result, $F$ is measured in units of eV/nm and $V$ in eV. The image charge factor is defined by $Q≡αfsℏc/4=0.36$ eV nm, where $αfs=1/137.036$ is the fine structure constant, and other terms have their usual meanings. A wave function of energy $E$ incident on an energy barrier $V(x)$ (in units of eV) then has an associated Gamow factor $θ[E(k)]$ and Kemble transmission probability defined by

in the simplest formulations used to analyze Fowler–Nordheim-type equations.⁵⁵ The “$o$” subscript on $θo$ reinforces that this is the conventional (unmodified) Gamow factor [and carries over to $Do(E)$]. In its place, a form $θ(E)=A+Bθo(E)$ described below may be substituted into $D(E)$ and “generalized” forms of each obtained. The relationship between energy $E$ and wave number $k$ shall be restricted to $E(k)=ℏ2k2/2m$⁠: thus, $θ$ can be expressed in terms of $k$ (the preferred dependence in the TMA methods given its relation to the definition of quantum mechanical current density $jk$ in Schrödinger’s equation³³) or $E$⁠, and therefore, below, the designations $θo[E(k)]$⁠, $θo(E)$⁠, and $θo(k)$ are interchangeable, with all referring to Eq. (1), depending on whether energy or momentum is under discussion. Forbes refers to $θo(E)$ as the “JWKB exponent” and sets the prefactor of $Pe−θo$ to unity,⁵⁶ or $P→1$⁠, in the definition of the Kemble form from the Fröman and Fröman⁵⁷ form $Pe−θo/(1+Pe−θo)$⁠. The Gamow factor accounts for sequential tunneling effects but not resonances,⁵⁸ which are inherently determined by the wave nature of the tunneling particle and, therefore, affected by abrupt transitions in the barrier. An example of the effects of abruptness is readily shown by using either the exact form of the Fowler–Nordheim triangular barrier $DFN(E)$ [Eq. (34) of Ref. 59, Eq. (2.19a) of Ref. 60, and Sec. 18.2 of Ref. 33] compared to the presence of an additional abrupt transition downstream as for an MIM barrier, or the inclusion of a defect well in the barrier itself^14,16,58 that introduces resonant tunneling (RT), as shown in Fig. 1: the narrowness of the well region and the distance of the right hand side drop-off would presumably not affect a Gamow factor approach to the evaluation of $D(E)$⁠, but clearly, the abrupt features associated with MIM and RT configurations introduce stark changes to the analytic FN $D(E)$ smooth curve. Properties of the SN barrier will cause the prefix $P$ in the Fröman and Fröman form⁵⁷ to depart from unity in general: as a first step, the generalized form $θ(E)=A+Bθo(E)$ is equivalent to replacing $P$ by a constant value.

To investigate the adequacy of the Kemble approximation for the SN barrier, Eq. (1) is compared to an exact TMA evaluation. The behavior of the exact analytic cases (principally, rectangular and triangular barriers⁴⁵) suggests that the reformulated Gamow factor, for which $Θ(k)≡ln⁡[(1/D(k))−1]$⁠, where $D(k)$ is the transmission probability numerically evaluated using TMA as per the conventions of Ref. 45, will bear an almost linear relationship to the simple Gamow factor away from the barrier maximum, as additionally suggested by the Fröman and Fröman form. The SN barrier for the TMA analysis is parameterized as

where $F=q|E|$ is related to the applied field $E$ and $Q/x$ is the image charge term. The maximum of the barrier occurs at $xo=Q/F$ and is $V(xo)=μ+Φ−4QF$⁠. Discretization of the barrier in Airy-TMA is shown in Fig. 2(a) for typical copperlike parameters (⁠$μ=7$ eV, $Φ=4.5$ eV, $F=6$ eV/nm, and $Q=0.36$ eV nm): unlike plane wave TMA, Airy-TMA requires few linear segments to characterize the barrier: although 19 points are used for visual smoothness, 15 points would not significantly degrade accuracy. The points are chosen such that their spacing in energy squared is equally distributed so as to favor more points near the apex; that is, $V(xj)=V(xN+1−j)=V(xo)[j/(N−1)]2$ for $j∈(1,…,N)$⁠, where $2N+3$ is the number of points defining the barrier. The two points just below the apex, however, are selected to better resolve the maximum of $V(x)$⁠. The associated TMA-determined Gamow factor $Θ(k)$ and transmission probability $D(k)$ are shown in Fig. 3 and compared to three approximations:

1. The baseline approximation $θo(E)$ of Eq. (1).

2. The refined approximation is defined by

   $θ(E)≡A+Bθo(E),$

   (3)

   where $A$ and $B$ are from a linear fit of the TMA $Θ(E)$ as a function of $θo(E)$ for $E<V(xo)$⁠.

3. The linear in energy approximation $θfn(E)≡θo(μ)−βF(μ)(μ−E)$⁠, where $βF(E)≡−∂Eθo(E)$⁠.

Approx. 1 is the standard JWKB approximation, Approx. 2 assumes that a linear relation exists between the baseline approximation and the TMA exact evaluation, and Approx. 3 is the approximation used by Murphy and Good in the development of the standard so-called Fowler and Nordheim equation for the Schottky–Nordheim barrier.^55,61,62

The roots $x±(E)$ of the SN barrier, defined by the relation $Vsn(x±)−E=0$ and Eq. (2), are

where $H(E)≡μ+Φ−E$ is the $x$-independent part of $Vsn(x)$⁠. The separation between the roots $L(E)=x+(E)−x−(E)$ is then

The roots $x±$ and $L$ both depend on the energy $E$⁠, but for simplicity in treatment, that dependence may not be shown but is presumed. The location of the maximum does not depend on $E$ and is

so that $Vsn(xo)=μ+Φ−4QF$⁠. Introduce the small parameter

and the ratio parameter

[observe that $y(μ)$ for the SN barrier is equivalent to $y=4QF/Φ$ from Murphy and Good⁶¹]. In terms of these parameters, the SN barrier is re-expressed as²² 

An integration over $x$ to obtain $θo(E)$ is numerically challenged by the behavior (shape) of the integrand, as numerical methods work best with smooth polynomial-like functions;⁶³ introduce $x=L[ϵ+sin2⁡s]$⁠, where $x−(E)=ϵL$⁠, and introduce the function $ho(s)$ defined by [compare to the more general $h(s)$ of Eq. (32)]

The integral of $ho(s)$ is well behaved: both it and its derivatives vanish at the integration limits. Therefore, evaluating the integral of $ho(s)$ using Simpson’s rule is highly accurate because of the Euler–MacLauren series method³³ even for a small number of summation terms, as discussed for Eq. (33). As a consequence, a rapid numerical method of evaluating $θo(E)$ is available that uses the simplest form of numerical integration (Simpson’s rule), is comparable in speed to the Forbes–Deane approximation, and is more accurate even for a small number of terms in the summation.

The SN Gamow function $θsn(E)$ is then expressible using the shape factor approach by

where $κ(E)$ is

Comparing the form^58,64

to Eq. (11) and using Eq. (5) for $L(E)$ allows the shape factor $σsn(y)$ to be quickly defined: it depends only on the dimensionless $y(E)$ of Eq. (8) and is

It is a weakly varying function with $y$⁠, with end points given by $σsn(1)=π/4$ and $σsn(0)=2/3$⁠. As a result, the shape factor method is very well suited to GTF methodology, as the primary variation of $βF(E)=−∂Eθo(E)$ is due to the dependence on $κ(E)$ and $L(E)$⁠, both of which are analytic, with $∂Eσsn(y)$ providing only a small modification.³⁰ To a good approximation,

where $z(y)≡(1−y)/(1+y)$⁠. The approximation is demonstrated in Fig. 4. Alternately, a small number of $σsn(y)$ values can be found exactly, and interpolation can be used.

Last, the $Q$ was given for a planar metallic surface, but the methodology above is indifferent to whether it is modified for semiconductors,⁶⁵ many body effects,^66,67 or curved surfaces:^35,44,58,68 the former will change $Q→(Ks−1)Q/(Ks+1)$⁠, and the latter two will change $Q/x$ to a more complex function that can be expanded to alter the effective work function and field in addition to contributing to the quadratic term considered below. As a result, $(Φ,F,Q)$ are treated as parameters that are easily understood in a planar metallic model, but which can be treated as effective parameters when curvature, many-body, and material properties demand and, therefore, which can make use of the methods here without further modification. The quadratic term, however, requires separate treatment.

The shape factor method^30,44,54 gives an alternate method to the Schottky–Nordheim functions for the treatment of the standard image charge barrier (alternately image-rounded barrier⁶⁹ or the SN barrier^70,71), but past developments using the shape factor approach have treated modifications to the barrier potential separately from the SN barrier. Here, a standardized theory is developed that has as an asymptotic limit the SN barrier, but that will be clearly generalizable to more complicated barriers such as MIM, depletion/accumulation, and interface barriers. For MIM barriers, in particular, a second image charge term must generally be included; if the insulator layer is sufficiently thin, then an infinite sum of image charge terms must be included, but as shown before, the shape factor method can be adapted to deal efficiently with that configuration.²² 

Let modifications to the SN barrier of Eq. (2) be of the form of a polynomial $Pn(x)$ of degree $n$⁠, or $Pn(x)=∑j=2nCjxj$ so that

where the subscript denotes the order of the polynomial $Pn(x)$ (and $Pn$ bears no relation to the Kemble factor $P$ above). Observe that although $j$ could start at 0 so that $C0=μ+Φ$ and $C1=−F$⁠, it is preferable not to do so and retain the special status of the linear terms explicitly, given the importance of the first two roots of the general barrier equation. The term containing $Q$ is the image charge contribution: in conventional units and for metals, $Q→Qo≡q2/16πϵ0=αℏc/4$ or 0.36 eV nm. It becomes $Q(Ks)=[(Ks−1)/(Ks+1)]Qo$⁠, where $Ks$ is the dielectric constant,⁶⁵ for semiconductors. By virtue of the tunneling problem, restriction may be made to polynomials such that $Vn(x)−E$⁠, where $E(k)$ is electron energy and $k$ is wave number, has only real roots $xj(E)$⁠: imaginary roots are associated with conditions for which the tunneling probability $D(E)≈1/{1+exp[θ(E)]}$ is driven to zero (infinitely thick barrier). The shape factor $σn(E)$ is then explicitly defined by

where $xo$ is the location of the closest local maximum to the surface and $(x−<x+)$ are the smallest two of the $xj$ roots defined by

In terms of the shape factor, the Gamow factor $θn(x)$ is then defined by

Observe that the form of $L(E)$ did not change, but that the expressions used for $x±(E)$ do. Thus, below, the $(±)$ subscripts on $x±$ shall be changed for the generalized barrier so as to designate roots by a subscript, or $xj$⁠.

It will be the convention below to retain $H$ as the height (in energy) of the barrier to the tunneling electron or $H(E)≡μ+Φ−E$ as in Sec. III: by virtue of the definition of $xo$⁠, it is seen that $xo$ is dependent only upon $F$⁠, $Q$⁠, and $Pj$⁠, but not on $μ$⁠, $Φ$⁠, and its dependence upon $E$ is only through the dependence of the coefficients $Pj$ in turn on $E$⁠. By the definition of polynomials with real roots, it is seen immediately that

In all that follows, Eq. (21) is the preferred form of the shape factor $σn(E)$ for further analysis. For numerical work, therefore, the task is to efficiently find the roots $xj(E)$⁠, the location of the turning points $xo$ and $xp$⁠, both defined by the solutions of Eq. (19) where the slope of the potential vanishes. The potential is assumed to be such that $Vn(x>xp)=Vn(xp)$ (constant).

For practical reasons, attention is now restricted to $n=2$ such that $P2(x)≡γx2$⁠: it is seen that $xV2(x)$ is a cubic: consequently, this is referred to as the “quadratic barrier model.” This is done for two reasons: first, the method to be developed below will be seen to be easily and straightforwardly extended to larger $n$ so that exposition is simplified by focusing on the $n=2$ case; and second, the roots $xj$ for $j>2$ will be large, making their contribution to Eq. (21) comparatively small. As a leading order approximation, therefore, quadratic $P2(x)≡γx2$ is adequate for the leading order approximation modified at most by a small factor to be given below: it is suggested to be, in general, sufficient for numerical work.

Cubic equations with three real roots are historically known as cases irreducibilis (irreducible case).⁷² The roots are sums of complex conjugates and so can be re-expressed in terms of trigonometric functions. For the quadratic potential barrier, the roots $xj$ and maximum/minimum locations $xo$ and $xp$ are specified by

where the root $xp$ follows the same equation as $xo$⁠, but other than observing that $xo$ and $xp$ lie between the first and second pair of roots, respectively, $xp$ is not important in the present discussion. The $xj$ roots satisfy the relations

and are directly given by

where $x2<x1<x3$⁠, an ordering dictated by the dependence of Eq. (24) on $j$ and not a consequence of preference. It is seen that $x−⇔x2$ and $x+⇔x1$⁠, and henceforth, the $(±)$ subscripts shall be discontinued for general barriers. The angles are defined by

and the factors $δ(E)$ and $y(E)$ are defined as

[compare Eq. (8)]. Both $δ$ and $y$ can be defined in terms of the roots $xj$ using Eq. (23),

where $⟨⋯⟩$ denotes average. The factor $y(E)$ is the ratio of the Schottky barrier lowering factor $4QF$ with the height $H$ of the barrier above the electron energy $E$ [Eq. (20) of Ref. 30]: for the special case of $E=μ$⁠, then, $y(μ)$ corresponds to the $y$ parameter of Murphy and Good.⁶¹ In computational parlance, Eq. (29) is vectorizable if the roots are specified: tables of lookup values, shown to significantly speed up execution when millions of emission elements are needed,³⁷ can, therefore, be economically calculated. For example, generation of a parabolic image charge barrier using $γ=0.225$ and $xj=[2.0,0.2,4.0]$ is shown in Fig. 5; also shown is the standard SN barrier (blue) for the same roots $x1$ and $x2$ but for a field of $F=Q/x1x2$⁠.

For $γ>0$⁠, then $H$ is constrained: for energies $E<V(xp)$⁠, the transmission probability is driven to zero because the barrier has become infinitely thick. Equations (25) and (26) can then be used to bound $E$⁠. For a given $H$⁠, the transmission probability vanishes when $γ$ exceeds a value $γmax$ for which $x1=x3$⁠, or, using Eq. (28), such that

For the purpose of analysis, though, specifying the third root $xj$ is not directly useful because the specification of three roots chosen independently will alter the strength of the image charge factor $Q$ and that is unphysical. It is, therefore, preferable to eliminate the largest root $x3$ in favor of $Q$ such that $x3≡Q/γx1x2$⁠: therefore, $σ$ will be constructed in terms of $x1,x2,$ and $Q$ for a given $γ$⁠. The SN barrier is recovered in the limit (⁠$γ→0$⁠) corresponding to $x3→∞$⁠. Ordered by size, the roots are $x2<x1<x3$⁠. Perform a change of variables as in Sec. III for the SN barrier such that $x→x2+Lsin2⁡s$⁠. It follows (suppressing the $n=2$ subscript on $σn$⁠):

where $Ho=V(xo)−E$ [the maximum of $V(x)−E$], $ϵ=x2/L$⁠, and $r≡L/x3$⁠. Observe that $ho(s)$ of Eq. (10) corresponds to $h(s)$ in the limit $r→0$⁠. All the roots and turning point locations have a dependence on $E$ because they depend on $H$⁠. Three limiting cases are analytic:

1. Triangular barrier: $Ks→1$ and $γ→0$ for which $x2→0$ and $ϵ→0$⁠. Then, $∫0π/2h(s)ds→4/3$⁠, $LQ/4Hox1x2→1/2$⁠, and therefore, $σ→2/3$⁠.

2. Parabolic barrier: $Ho→0$ (barrier becomes a parabola), $x1≈xo+L/2$⁠, and $x2≈xo−L/2$ as $L→0$⁠. $∫0π/2h(s)ds→(π/4)L(x3−xo)/xox3$ and $LQ/4Hox1x2→xox3/[L(x3−xo)]$⁠, and therefore, $σ→π/4$⁠.

3. Depletion barrier: $Ks→1$⁠, $Ho=V(xo)−V(xp)$⁠, then $x3=x1$⁠, $x2→0$⁠, and $xo→0$ (concave up parabolic barrier without image charge). $LQ/4Hox1x2→1/2$⁠, and therefore, $σ→1/2$⁠.

The form of Eq. (31) is advantageous because the integral involving $h(s)$ varies over a narrow range, and the remaining terms are obtained through direct function evaluations. Because both $h(s)$ and its derivatives vanish at the integration limits, an Euler–MacLauren approximation makes Simpson’s rule a good approximation.³³ Consequently, the summation

where $sj=jπ/2(N+1)$⁠, is numerically expedient because it can use vectorized algorithms, making the evaluation of $σ$ for small $N$ computationally rapid. The behavior of $h(s)$ for a constant $Q=0.36$ eV nm and fixed roots $(x1,x2)=(2,0.2)$ in nm, but varying $x3$⁠, which serves to change all $(γ,F,H)$⁠, is shown in Fig. 6.

The evaluation of Eq. (31) exactly using full numerical integration is compared to using Eq. (33) for various $N$'s in Fig. 7 for three cases characterized by the conventional representation with $Vo=μ+Φ$ and $Q=(Ks−1)Qo/(Ks+1)$ to describe the barriers. The barriers are characterized by the four conventional parameters $(Vo,F,Q,γ)$ and turning points $(xo,xp)$ in Table I. The associated barriers are shown in Fig. 7(a) and the corresponding $σ(E)$ for $0≤E≤V(xo)$ in Fig. 7(b). In the latter, symbols denote the exact evaluation using Eq. (31), and the lines (dashed, thin, thick) denote the $N$ in increasing order in the approximation Eq. (33). Smooth, image charge rounded barriers allow for even smaller $N$ resulting in good approximations, but even a difficult depletion barrier with only a minuscule image charge modification (Case 1) is well approximated for a manageable size of $N$⁠. As a result, the Gamow $θ(E)$ factors can be quickly and accurately evaluated without resorting to numerical integration.

To evaluate current density using a GTF-based approach, $θ(E)$ is best represented by a cubic in energy. Doing so allows $βF(E)$ to be found directly using standard polynomial methods for derivatives. Moreover, the location of the current density integrand is then rapidly found using optimized search algorithms, allowing for the reformulated GTF method³⁰ to be directly used in evaluating the GTF function $N(n,s)$⁠. This line of inquiry will be the subject of a separate study. In the remainder of the present work, representative temperatures for nanotubes and fibers drawing significant current are investigated: from those temperatures and the behavior of the generalized shape factors and Gamow factors, a preliminary assessment of the changes incurred by quadratic barriers, due to space charge or curvature, can be evaluated. Here, it is observed that the primary consequences of general quadratic barriers will be to marginally reduce the thermal contributions due to small modifications to $V(xo)$⁠. The impact on field emission contributions will be proportionately greater due to the enlargement of the tunneling barrier and, therefore, a suppression of the tunneling current. As a result, the onset of the thermal-field regime,³⁰ where $n=βT/βF≈1$⁠, is expected to be moved to higher temperatures. A quantification of that assessment shall also be the subject of a separate study.

In a prior study, the temperature variation along a damaged emitter was shown to increase.²⁵ When conductors undergo constrictions, either intentionally as a result of fabrication¹⁸ or due to operational processes,^25,73 an increase in the temperature due to a smaller cross section is anticipated^74,75 as cylindrical emitters taper to a “point” allowing it to be modeled as a conical emitter. The temperature profile of a cone shaped emitter is, therefore, expected to depart from a cylindrical emitter along the length of the cone using physics similar to the cylindrical study.²⁵ When large currents were drawn from carbon fiber emitters,⁷³ loss of material caused the wire diameter to thin and the apex to taper and give evidence to the large temperature excursion hypothesis. To model the conductivity, one must treat the change in a cross-sectional area as a gradually decreasing process.

The images of Fig. 8 provide a model of the cross-sectional area at the tip of the cone, where the depth and extent of perforation based on the frequency of visible holes and color contrast in the image suggest a plausible cone model of the damaged wire with a certain half-angle $θ$⁠. A histogram of the pixel data used to create the gray scale contour map quantifies the lightest parts of the image and is shown in Fig. 9. A summation over the lightest parts of the histogram is then taken to be proportional to the cross-sectional area of the wire at that point in the conical model. As evident in Fig. 10, the perforations become more extensive near the apex. The depth and frequency of perforation at varying points along the fiber were assessed based on the frequency of visible holes and color contrast of images from the prior study (see Fig. 1 in Ref. 25). Using an image of the tip of a fiber shown in Fig. 8, the variation in lightness of the pixels on a gray scale can be used to estimate the cross section such that a sum of the lightest pixels is proportional to the cross-sectional area of the fiber at that location. Estimates were also based on measurements of the depth and perforation frequency along the wire from associated images. The cross-sectional area estimates allowed for the modeling of damage as a cone with a half-angle of approximately $20°$ as shown in Fig. 10. As the frequency and magnitude of the holes increased, much of the exterior part of the damaged tip was omitted due to the capacitorlike character of ledges and openings in the geometry. Omission of the exterior perforations allowed for the interpretation of the current carrying portions of the wire as a solid conical core. The relationship between the temperature and the cone angle was investigated by performing simulations for the different cone half-angles of Fig. 11. The resulting temperature profiles are shown in Figs. 12 and 13.

Observe that these evaluations were performed to assess heating at the apex of a long, cylindrical carbon fiber for which the apex was undergoing heating and deformation and likely sublimation. The same conical evaluations, however, can be performed for actual conical emitters, such as the 12CaO-7Al$2$O$3$⁠:4e$−$⁠, or C12A7, electride emitters characterized by a wider half-angle conical shape closer to the blue case of Fig. 11 (e.g., Fig. 1 in Ref. 18) but possessing small nanometer-scale “cages” and is under investigation as a hollow cathode for spacecraft plasma propulsion.^18,76 Such emitters are possibly characterized by conditions closer to the “copper” parameters herein if the conductivity is greater.

The evaluation of heating for a conical emitter must generalize the one-dimensional differential equation modeling temperature variation along a cylinder^25,77,78 by extending the diffusive term (dependent on $κ∂x2T$ when the thermal conductivity $κ$ is independent of temperature, as it is in the Drude model¹⁸) to radial coordinates. Doing so significantly increases the complexity of the temperature evaluation for when current is compressed along a conical conductor. Note, however, the study of Kyritsakis et al.,⁷⁹ who found in favor of the applicability of 1D models even to 3D structures experiencing thermal runaway. Although heat diffusion through the base of the conical emitter cannot, in general, be neglected in such evaluations, we observe the following. First, a simple analytical model of the temperature along a cylindrical conducting wire based on solutions to a variant of the harmonic oscillator equation and which neglects the radiative cooling term already predicts substantial temperature rises near the apex.²⁵ That simple model can be amended to include radiative heat loss, and that study demonstrates even large temperature excursions. Here, temperatures in the cone will still be higher due to increased current constriction. Second, in addition to the study by Kyritsakis et al., a detailed heat study of conical field emitters on both flat conducting surfaces as well as tip-on-post geometries was conducted by Ancona⁸⁰ (who also included Nottingham heating near the emission site as did Kyritsakis et al., heating that can transition to Henderson cooling at elevated temperatures⁷⁸) who found that even for metallic (molybdenum) emitters, small cone angles were associated with a relatively poor heat sinking capability to the base, an effect exasperated significantly when the conical emitter was atop a post—a geometry more in line with the present study. Third and finally, Kyritsakis et al. found sufficient thermal energy present to dynamically elongate the emitter during high current/high temperature operation that contributed to thermal runaway. Present purposes seek an indicative, not an exact, account and, therefore, investigate instead a representative model for which heat losses through the base of the conical emitter are neglected so as to establish an upper limit to the temperature excursions that should be anticipated, even if those excursions are higher than expected values if heat diffusion were properly accounted for, with the separation between actual and upper limit expected to increase for good thermal conductors and large cone angles.

Temperature profiles for two types of conical emitters were evaluated, one taken to be carbonlike and the other to be metallic (copper) for comparison. Temperature along the length of the cone was evaluated assuming steady state current flow so that no net heat was supplied to any element of the cone. For the region shown in Fig. 11, balancing only Joule heating and radiative cooling processes (that is, neglecting the diffusive heat loss term through the base of the cone) for a cone with half-angle $θ$ given by

where $z$ is the vertical distance along the cone (with $z=0$ at the apex); $Ac=πr2$ is the average cross-sectional area of a segment; $dAs=2πrdl$ is the surface area of circumference $2πr$ and thickness $dl=dr2+dz2=dz/cos⁡θ$⁠, where $tan⁡θ=ro/h$ is the cone angle; $ro$ is the base radius, and $h$ is the cone height. Electrical parameters are $ρ$ is the resistivity of the material, $I$ is the current through the fiber and the same for each segment, $ϵ$ is the emissivity, $σB$ is the Stefan–Boltzmann constant, and $T$ and $To$ are the temperature of the segment and environmental background temperature, respectively. Including diffusive heat loss through the inclusion of a term going as $κ∂x2T$ lowers the value of $T$ satisfying Eq. (34), but at the cost of the rapid analytical treatment pursued here; its inclusion is to be taken up separately. Manipulating gives

Resistivity as a function of temperature for metals and metal-like materials is taken to be

(observe that here, the sign of $α$ changes depending on a metal or a semiconductor), which in combination with Eq. (35) gives

Temperature roots of this equation then determine the temperature profiles shown in Figs. 12 and 13 and discussed below. For each radii value $r$⁠, there is an associated vertical distance $z=r/tan⁡θ$ along the cone. This procedure is indicative but is undermined near the apex, where current is increasingly pushed through a vanishing cross section: in actuality, processes that ablate and blunt the tip are expected to intercede before that.⁷⁹ 

The resistivity of a cone for three different half-angles was considered. The $20°$ angle is representative, but the $30°$ and $10°$ were investigated as well to assess sensitivity. As the cone half-angle decreases, the temperature near the apex is significantly higher than for larger half-angles for both carbon and copper cases: the cross-sectional area decreases more gradually in a cone with a smaller half-angle and affects the diffusion of heat.

As seen in Figs. 12 and 13, the upper limit temperature of a conical emitter for the specific analytical model of Eq. (37) steeply increases in proximity to the cone tip, as anticipated, for a steady state current taken to be 10 mA, a representative value for these fiber emitters. The variation of the temperature strongly depends on the material parameters of Table II. Smaller values of $α$ are indicative of more conductive materials. For carbon, few configurations allow a temperature less than 1000 K, while for copper, temperatures exceeding 1000 K only occur within 5–20 $μ$m of the tip due to the faster heat conduction associated with copper behavior. As the sublimation temperature of carbon is high, most of the conical emitters resist sublimation down to radii below 20 $μ$m. The melting temperature of a metal is much lower (e.g., for copper, it is 1358 K). Nevertheless, at 10 $μ$m radii, most of the fiber is expected to ablate material and become blunted owing to the high-field enhancement of a fiber and the high current density from the emission sites, as seen in Fig. 8 where damage to the carbon fiber is evident.

In our prior study,²⁵ both the numerical and analytical models predicted large temperatures beyond $To$ by the location of the base of the conical region considered here. Therefore, the overoptimistic declines in Figs. 12 and 13 in the direction of room temperature values are a measure of the importance of solving the full temperature gradient equation, including the diffusive term $κ∂x2T$⁠. The intent here, however, is to show the sensitivity of the temperature along the cone due to the cone angle by way of a simple model sufficient for specifying upper bounds, where the variation in the constants listed in Table II can have a massive affect on temperature. If the emissivity listed for copper is $ϵ=1.0$ instead of $ϵ=0.1$⁠, the temperature nearly doubles. The values given in Table II are arbitrary values. Figures 12 and 13 are dependent on the values listed in Table II; therefore, the temperature of the conical emitter at different lengths along the wire can vary. If the upper limit temperature is comparable to the actual temperature when the diffusive term in the thermal gradient equation is retained, then the relative variation in the temperature shown here suggests a reevaluation of the probability of sublimation or melting of a conical emitter, a finding supportive of and supported by Refs. 79 and 80. Under different conditions, the functions in Figs. 12 and 13 can be recalculated with variables that are representative of the material in question giving a greater feel for the stability of the emitter as it heats. Therefore, it is reasonable to use the equations provided in this section to model the temperature of the cone under the constraints that they are used with the appropriate parameters based on the material at hand and under limited conditions where the thermal diffusion through the base is constrained.

At very high temperatures, thermal-field emission occurs, causing the presumed constancy of temperature along the cone to be progressively undermined, although field emission can dominate (the neglected Nottingham heating contributions that may be significant at lower temperature can turn to Henderson cooling at a higher temperature). To estimate this effect, the total current from differential surface elements must be estimated using the revised GTF equations accounting for the barrier modifications associated with the curvature considered in Sec. IV. Furthermore, depending on the desired affect, the half-angle developed and material parameters for the conical emitter can be used to alleviate or create higher temperatures depending on their use.

Both nanotube and conical field emission sources can run at elevated temperatures, making a description of their thermal-field behavior complicated. Further complications are gradients in the temperature along the emitter, separate regions of the emitter exhibiting thermal, field, or thermal-field behavior, variations along the emitter processes, such as space charge or curvature, which alter the barrier shape and, therefore, change the barrier height $μ+Φ$ and the tunneling width $L(E)$ of it, and material processes that can change the nature of the image charge contribution itself. Many of these contribute simultaneously and complicate the use of fast emission models in PIC and MD codes. In the present study, we have sought to

• show how exact treatments can lead to the development of generalized Gamow factors,

• develop computationally rapid means of obtaining the Gamow factor for variations,

• show how a generalized Gamow factor can be used in the GTF equation for thermal-field conditions, and

• assess how a tapered structure can contribute to the substantial elevation of the presumed temperature of the emission cone.

Processes that reduce $Q$ or increase $L(E)$ serve to diminish field emission components; processes that increase temperature increase the size of the TF region on the emitter and alter the distribution (and velocity) of launched electrons off the surface. Such features are difficult in the modeling of vacuum nanoelectronic devices or thermal-field emitters, which are intentionally operated at elevated temperatures. In the present study, developments designed to enable flexible and computationally rapid emission models to provide launch conditions in advanced beam codes have been treated, and indications for how they can be used or extended were given.

K.L.J., A.S., and M.O. gratefully acknowledge support from the NRL Naval Innovative Science and Engineering (NISE) program.

The data that support the findings of this study are available within the article.

The standard Schottky–Nordheim barrier results are obtained in the limit $γ→0$⁠, but the process requires some care. The small $γ$ limits of the roots $(x1,x2)$ are most easily handled by considering the combinations $x1±x2$ and $xo$⁠. First,

To leading order, $φ=δ[(3/4)(1−y2)]1/2$⁠, and therefore,

Second,

To leading order in $φ$ again, or $cos(φ/3)→1$⁠,

Last, for $xo$⁠, use $χ≈2γ3Q/F3=yδ/3$ to find in the limit of small $γ$⁠,

As $δ=3γH/F2→0$ in Eqs. (A2), (A4), and (A5), the Schottky–Nordheim, or the standard image charge barrier, results [namely, Eq. (6) and $xp→∞$] are recovered.