Field emitter electrostatics: a review with special emphasis on modern high-precision finite-element modelling

The electrostatics of field electron emission (FE) is most commonly (and most easily) discussed by using classical-conductor models in which the emitter and relevant parts of the enclosing experimental system are represented by smooth geometrical surfaces (the 'conductor surfaces') on which the classical electrostatic potential ${\varPhi}$ or (in some contexts) the conventional classical electrostatic field $\boldsymbol{E} \; \left[{=}{-}\mathbf{grad}{\,}{\varPhi}\right]$ is specified. Inside the conductor surfaces (i.e. on the conductor side) the field is taken as zero; outside, the classical electrostatic potential and field are taken to obey Maxwell's equations in the form of Laplace's equation. Classical-conductor models work best for metals, but are commonly applied to FE from non-metals as a first approximation.

We now prefer to call E the 'electrostatic field' rather than the 'electric field', because the latter term can be understood to also cover the optical E-fields generated by lasers. Lasers are increasingly used in FE contexts. However, as indicated above, applied electrostatic fields and applied optical E-fields can generate different physical effects. In particular, the processes of field enhancement are physically different in the two cases.

In most FE situations, it is customary to disregard work-function variations across the emitter and system surfaces (and the resulting patch fields). Thus, in the absence of current flow along an emitter, the potential ${\varPhi}$ in the part of the conductor surface that defines an emitter is treated as constant. Further, if the surfaces of all system electrodes are treated as having the same work function as the emitter, as is usually done for simplicity, then the classical potential difference $\Delta {\varPhi}$ between a counter-electrode and the emitter electrode is equal to the voltage V applied between them. In this case, the technical distinction between $\Delta {\varPhi}$ and V, which only rarely has any significant quantitative consequences, can be disregarded when formulating basic theory.

The above remarks apply to classical-conductor models. However, with both field ion and field electron emission there are experimentally observable atomic-level phenomena. There exist both an older body of atomic-level 'electrostatic' theory (developed in conjunction with the theory of how the field ion microscope works), and an increasing modern body of atomic-level electrostatic theory–mainly in the context of FE (and particularly in the context of FE from CNTs), and mainly implemented via density functional theory (DFT). Deep issues arise about how to properly connect classical-conductor and modern atomic-level theories of electrostatics. These issues, and atomic-level theories of electrostatics in general, are outside the scope of this review, which focuses on the accuracy of classical-conductor electrostatic models. However, the material in this review should assist in later discussion of atomic-level electrostatic models.

All quantities and equations in this article are defined using the modern international system of equations that has the vacuum electric permittivity ε₀ in Coulomb's Law. Since 2009, this system has been called the 'International System of Quantities' (ISQ). Thus, voltages and classical electrostatic potentials are measured in volts, and fields are measured in V m⁻¹ (or a related unit with the same dimensions, such as V nm⁻¹). Most published FE articles now use the ISQ, though some older FE textbooks use Gaussian system equations to formulate aspects of FE theory.

In field emitter physics, it is customary to take the direction away from the emitter into the vacuum as the positive direction. This convention implies that the electrostatic fields that induce field electron emission are negative in sign. However, the 'fields' that go into the related equations for local emission current density (LECD) (J_L) and current (I) need, by convention, to be positive in sign. Customary practice in FE literature is to use the symbol E or F to denote a positive quantity that is the absolute magnitude of electrostatic field. We think it clearer here to denote this positive quantity by $|E|$. However, we use the customary convention that the symbols J and I denote positive quantities that are the absolute magnitudes of current density and current, respectively.

In FE theory, the main role of an electrostatic field is as a parameter in an emission equation. A convenient starting point is the zero-temperature version of Murphy–Good (MG) FE theory [15], which was developed in 1956 in order to correct serious errors in the original 1928 theory of Fowler and Nordheim and/or in directly related 1928 and 1929 papers. (For a modern derivation of MG FE theory, see [16].)

Notwithstanding this, many modern FE technological papers use a simplified version of the original 1928/29 FN FE theory, rather than MG FE theory. It has been argued [17, 18] that this is weak scientific practice, and that using MG FE theory would lead to better science and to better value for research funders.

Both FN FE theory and MG FE theory are based on a simple underlying physical model, originally introduced in the 1920s, that (a) disregards the existence of atomic structure and the role of atomic-level wave-functions on the emission process, and (b) models a field electron emitter as a Sommerfeld-type free-electron conductor with a smooth planar surface of large lateral extent. This has been called smooth planar metal-like emitter (SPME) methodology.

In MG theory, FE takes place by tunnelling from emitter electron states near the emitter Fermi level, through a barrier modelled as a planar image-rounded barrier–often now called a Schottky–Nordheim (SN) barrier [16]. In the zero-temperature theory, the resulting LECD $J_{\textrm{L}}^{\,\textrm{MG0}}$ is a function of the local work function φ and of the absolute magnitude $|E_{\textrm{L}}|$ of the constant electrostatic field above the planar emitter surface.

For present purposes it is convenient to write the zero-temperature 1956 MG FE equation in the 'linked' form

$$\begin{equation} J_{\textrm{L}}^{\,\textrm{MG0}} = {{{\textrm{t}}}_{\textrm{F}}^{-2}} J_{\textrm{kL}}^{\,\textrm{SN}}, \end{equation} \tag{ 1 }$$

$$\begin{equation} J_{\textrm{kL}}^{\,\textrm{SN}} \equiv a {\phi}^{-1}|E_{\textrm{L}}|^{2} {\exp{[-{{\textrm{v}}}_{\textrm{F}}} b {\phi}^{3/2}/|E_{\textrm{L}}|]}, \end{equation} \tag{ 2 }$$

where a and b are the Fowler–Nordheim constants [19], and ${{\textrm{v}}}_{\textrm{F}}$ and ${{\textrm{t}}}_{\textrm{F}}$ are particular values (appropriate to a SN barrier defined by φ and $|E_{\textrm{L}}|$) of the field emission special mathematical functions ${{{\textrm{v}}}}(x)$ and ${{{\textrm{t}}}}(x)$ [5], where x is the Gauss variable, i.e. the independent variable in the Euler/Gauss Hypergeometric Differential Equation. The function $J_{\textrm{kL}}^{\,\textrm{SN}}$ is defined by equation (2) and is called the local kernel current density for the SN barrier.

Obviously, many real emitters do not have large flat emitting surfaces, but have the shape of a pointed needle or a rounded post, or are otherwise 'pointy'. We call these pointform emitters. For pointform emitters that are 'not too sharp', MG FE theory can be applied by using the so-called 'planar emission approximation'.

Consider a location 'L' in the emitter surface, and assume the kernel current density at 'L' is given by equation (2). Integration over the emitter surface yields the notional emission current (for the SN barrier)

$$\begin{equation} I_{\textrm{n}}^{\,\textrm{SN}} = \int_{\Sigma} J_{\textrm{kL}}^{\,\textrm{SN}} \textrm{d} \Sigma .\end{equation} \tag{ 3 }$$

$I_{\textrm{n}}^{\,\textrm{SN}}$ also depends on the emitter shape and the assumed common work-function value, but these dependences are not explicitly shown.

The planar emission approximation is poor if the emitter apex radius is 'too small'. This is because the correct curved-surface expressions for the barrier transmission probability and the electron supply then diverge significantly from those applicable to planar geometry. 'Too small' is often taken to mean 'radius less than 10–20 nm'—though recent work [20] suggests that 'less than 100 nm', or even 'less than 1 µm', might be a better criterion, depending on what level of LECD precision is required (the 1 µm criterion corresponds to 5%). Nevertheless, the planar emission approximation has been very widely used in practice (possibly often without recognition).

After this integration, it is convenient to choose a characteristic location 'C' on the emitter surface where the local field magnitude (and hence the LECD) are particularly high, and then define a notional emission area, using 'C' and the SN barrier, $(A_{\textrm{nC}}^{\,\textrm{SN}})$ by

$$\begin{equation} A_{\textrm{nC}}^{\,\textrm{SN}} \equiv I_{\textrm{n}}^{\,\textrm{SN}}/J_{\textrm{kC}}^{\,\textrm{SN}}. \end{equation} \tag{ 4 }$$

In modelling, where cylindrical symmetry and uniform work function are nearly always assumed, 'C' is usually taken at the emitter apex.

To the best of our knowledge, at no stage in the last 90 years or so since the 1928 Fowler and Nordheim paper [21] have there been any comparisons of theoretical and experimental LECD values that were reliable, precise and decisive. Further, it is clear that theories that disregard the existence of atoms (which SPME methodology does) are not physically complete. Hence, it can be argued that FE theory is still in a pre-scientific state, and that it must be assumed that $I_{\textrm{n}}^{\,\textrm{SN}}$ is not an accurate estimate of experimentally measured current—even for so-called 'electronically ideal' FE systems where the measured current is determined by the emission physics and the system geometry alone, and is not influenced by any so-called 'system complications' (see section 2.4).

To deal formally with this situation, for electronically ideal systems, a prediction uncertainty factor (previously called a 'knowledge uncertainty factor') has been introduced (e.g. see [5]). The treatment here aligns with that in [18], is more complete than given previously, and supersedes earlier discussions.

We denote the real physical LECD at a given surface location 'L' by $J_{\textrm{real,L}}$ and take it as given by an equation of the form

$$\begin{equation} J_{\textrm{real,L}} = {\lambda_{\textrm{kL}}^{\,\textrm{SN}}} J_{\textrm{kL}}^{\,\textrm{SN}}, \end{equation} \tag{ 5 }$$

where the LECD prediction uncertainty factor $\lambda_{\textrm{kL}}^{\,\textrm{SN}}$ is defined by this equation. For simplicity, the pre-factor ${\textrm{t}}_{\textrm{F}}^{-2}$ in equation (1), and the temperature-dependent correction factor that appears in the finite-temperature version of MG FE theory, are 'swept up' into $\lambda_{\textrm{kL}}^{\,\textrm{SN}}$.

Values and functional dependences of $\lambda_{\textrm{kL}}^{\,\textrm{SN}}$ are not known. At present, the best guess of Forbes [22] is that it is a function of field and lies in the range $0.005 \lt \lambda_{\textrm{kL}}^{\,\textrm{SN}} \lt 14$.

When modified in this way, MG theory has been called [23] Extended Murphy–Good (EMG) FE theory.

For an emitter of given 'configuration' (i.e. of given shape and of given work-function distribution over its surface Σ), a true model emission current $I_{\textrm{tm}}^{\,\textrm{SN}}$ would be given by integration, and could be put in the form shown:

$$\begin{equation} {I_{\textrm{tm}}^{\,\textrm{SN}} = \int_{\Sigma} J_{\textrm{real,L}} \textrm{d} \mathit{\varSigma}} = \lambda_J I_{\textrm{n}}^{\,\textrm{SN}} = \lambda_J A_{\textrm{nC}}^{\,\textrm{SN}}{\;}J_{\textrm{kC}}^{\,\textrm{SN}}. \end{equation} \tag{ 6 }$$

Here, λ_J is an LECD–related prediction uncertainty factor of the same general kind as $\lambda_{\textrm{kL}}^{\,\textrm{SN}}$. λ_J will depend on many things, and we do not think it helpful here to introduce complicated notation about this: λ_J is basically a 'place-holder' indicating that uncertainty exists.

We now take into account that the true configuration of a real emitter may be different that assumed during modelling. This introduces a second uncertainty factor (λ_EM), called here the emitter-model prediction uncertainty factor, and results in a formal formula for the predicted emission current $I_{\textrm{p}}^{\,\textrm{SN}}$:

$$\begin{equation} {I_{\textrm{p}}^{\,\textrm{SN}} = \lambda_{\textrm{EM}} I_{\textrm{tm}}^{\,\textrm{SN}} = \lambda_{\textrm{EM}} \lambda_J A_{\textrm{nC}}^{\,\textrm{SN}} J_{\textrm{kC}}^{\,\textrm{SN}}}. \end{equation} \tag{ 7 }$$

If a leakage current exists, in a FE system that is otherwise electronically ideal, this can be treated as equivalent to a small change in λ_EM.

This formula can be simplified by defining a new parameter, $A_{\textrm{fC}}^{\,\textrm{SN}}$, the characteristic formal emission area as defined using 'C' and the SN barrier, by

$$\begin{equation} A_{\textrm{fC}}^{\,\textrm{SN}} = \lambda_{\textrm{EM}} \lambda_J A_{\textrm{nC}}^{\,\textrm{SN}}. \end{equation} \tag{ 8 }$$

Equation (7) can then be written

$$\begin{equation} I_{\textrm{p}}^{\,\textrm{SN}} = A_{\textrm{fC}}^{\,\textrm{SN}} J_{\textrm{kC}}^{\,\textrm{SN}} . \end{equation} \tag{ 9 }$$

We now assert that, in EMG FE theory, for an electronically ideal emitter, a formula for the actual experimentally measured current $I_{\textrm{m}}^{\textrm{EMG}}$ is

$$\begin{equation} I_{\textrm{m}}^{\textrm{EMG}} = A_{\textrm{fC}}^{\,\textrm{SN}} J_{\textrm{kC}}^{\,\textrm{SN}} = A_{\textrm{fC}}^{\,\textrm{SN}} a {\phi}^{-1}|E_{\textrm{C}}|^{2} {\exp{[-{\textrm{v}}}_{\textrm{F}} b {\phi}^{3/2}/|E_{\textrm{C}}|]} . \end{equation} \tag{ 10 }$$

In this equation, both the measured current and the kernel LECD are well defined in principle, so if (a) the work-function value φ is well known, and (b) the relationship between measured voltage and characteristic local field is well known (or can be determined experimentally), then an 'experimental' value for formal emission area can be extracted from experimental current–voltage measurements.

Note that, for a given emitter of presumed known configuration, the notional emission area is the area that would be predicted by theory, but the formal emission area is the parameter extracted from experiments. As things stand at present, even for a hypothetical emitter of known well-defined configuration, it would not be possible for theoreticians to accurately predict the 'area' value that would be measured by experimentalists, and it would not be possible for experimentalists to accurately measure an 'area' value predicted by theoreticians.

To ameliorate this problem and put FE theory onto a proper scientific basis, one needs much better knowledge (than currently exists) of the behaviour of all four of $A_{\textrm{fC}}^{\,\textrm{SN}}$, $\lambda_{\textrm{kL}}^{\,\textrm{SN}}$, $A_{\textrm{nC}}^{\,\textrm{SN}}$ and (ideally) λ_EM. In particular, it would be helpful to know how the notional emission area depends on the assumed work-function value and on emitter shape, particularly on the apex radius, and (using a slightly modified version of the theory) on temperature. An important new role for electrostatic simulations using finite element methods is to help explore questions of this kind. Some useful progress has been made [4, 24–26], but much more remains to be done.

A FE system is defined to include all aspects of the experimental system that can affect the measured current–voltage [$I_{\textrm{m}}(V_{\textrm{m}})$] relationship, including: (a) emitter composition, geometry and surface condition; (b) the mechanical, geometrical and electrical arrangements in the vacuum system; (c) all aspects of the electronic circuitry and all electronic measurement instruments; (d) the emission physics; and (e) all relevant physical processes that might be taking place (for example, the generation of field emitted vacuum space-charge, Maxwell-stress-induced reversible changes in emitter geometry, and adsorbate-atom dynamics).

The measured voltage V_m (normally measured at the high voltage generator (HVG)) is defined as the voltage applied between the counter-electrode side of the system and the emitter side of the system. Thus, in field electron emission the voltage V_m is positive in sign.

In general, the relationship between the electrostatic field E_L at location 'L' in the emitter's electrical surface and the measured voltage V_m can be written in the form

$$\begin{equation} E_{\textrm{L}} = -V_{\textrm{m}}/\zeta_{\textrm{mL}}, \end{equation} \tag{ 11 }$$

where ζ_mL is the 'measured' local voltage conversion length (LVCL). It is usually convenient to give LVCLs in nm.

For the analysis of measured $I_{\textrm{m}}(V_{\textrm{m}})$ characteristics, the relationship between V_m and the emitter apex field E_a is of particular interest, since this can be combined with equation (10) (with 'C' = 'a') to yield a formal expression for the $I_{\textrm{m}}(V_{\textrm{m}})$ characteristic.

An FE system is described as electronically ideal if, for all surface locations 'L', the parameter ζ_mL is constant, independent of the measured current or voltage, or of time. For ideal systems, the symbol ζ_mL can be replaced by ζ_L, and this 'ideal LCVL' is simply called the 'local VCL'. (Earlier treatments did not make this distinction, and the symbol ζ_L was used for both meanings.)

For these ideal systems, it is better to replace equation (11) by the 'ideal' or 'purely electrostatic' equation

$$\begin{equation} E_{\textrm{L}} = -V_{\textrm{M}}/\zeta_{\textrm{L}}, \end{equation} \tag{ 12 }$$

where V_M is the voltage between the counter-electrode and the emitter electrode; V_M has sometimes been called the 'macroscopic voltage' but is called here the inter-electrode voltage.

Strictly, electrostatic discussions should be formulated in terms of fields and electrostatic potentials, but it is often more convenient—and closer to actual practice in the literature—to formulate discussions as a relationship between fields and voltages, and ignore the distinction between voltage and classical electrostatic potential difference. (If the simplifying assumption of uniform local work function is made, there is no numerical distinction, anyway.)

In some real systems the inter-electrode voltage V_M is not equal to the measured voltage V_m, because there is significant series resistance in the current path between the HVG and the front (vacuum facing) surface of the emitter electrode. In particular, this can happen if an emitting protrusion stands on a poorly conducting semiconductor substrate. Thus, it is essential to make the distinction between ζ_mL and ζ_L.

More generally, many real FE systems are not electronically ideal, due to the occurrence of one or more 'system complications'. Possible complications include, amongst others: (a) significant series resistance in the conducting path between the protrusion and the HVG, as just discussed; (b) work-function changes, due to current-related heating (Joule and/or Nottingham heating) and resulting desorption of surface adsorbates; (c) geometrical emitter-shape changes due to Maxwell stress; (d) shape changes due to electroformation or emitter erosion; (e) effects due to field emitted vacuum space-charge (FEVSC); (f) current dependence in LVCLs, and in the related field enhancement factors defined below, due to voltage loss along the emitter; and (g) (for semiconductors) field penetration into the emitter surface.

A practical validity check, called an orthodoxy test [27], has been created that can indicate whether or not it is very probable that a given set of measured current–voltage characteristics relate to a FE system that is electronically ideal.

In FE literature, the apex ('a') version of equation (12) has been written in several equivalent forms, namely (using the notation preferred here)

$$\begin{equation} (-E_{\textrm{a}}) = \beta_{\textrm{a}} V_{\textrm{M}} = V_{\textrm{M}}/{\zeta_{\textrm{a}}} = V_{\textrm{M}}/{k_{\textrm{a}} r_{\textrm{a}}}. \end{equation} \tag{ 13 }$$

The parameter β_a is defined as $1/\zeta_{\textrm{a}}$, and is also written 'β' and 'β_V ' in the literature. These 'betas' are conveniently measured in m⁻¹. The symbol 'β' has been widely used in older literature (e.g. [28]) but clashes with the widespread modern use of β to denote the dimensionless field enhancement factor denoted below by γ.

Confusingly, the parameter β_a in the equation $|E_{\textrm{a}}| = \beta_{\textrm{a}} V_{\textrm{m}}$ is also sometimes called a 'field enhancement factor'. This dual use of the symbol β and the term 'field enhancement factor' can and sometimes does give rise to misunderstandings; hence, we now prefer to use ζ_a. More generally, our view is that a parameter measured in nm is easier to interpret than its reciprocal.

When equation (11) is rewritten in the form $V_{\textrm{m}} = \zeta_{\textrm{ma}} |E_{\textrm{a}}|$, it is readily seen that ζ_ma is a measure of how easy it is to 'turn an emitter on' in a given FE system. If, say, the emission onset value of the field magnitude $|E_{\textrm{a}}|$ is 2 V nm⁻¹, then a ζ_ma-value of 10 nm implies a turn-on voltage of 20 V, whereas a ζ_ma-value of 100 nm implies a turn-on voltage of 200 V, etc. Note that ζ_ma and ζ_a are characterization parameters, not physical distances.

In equation (13), r_a denotes the emitter apex radius, and the dimensionless parameter $k_{\textrm{a}} = \zeta_{\textrm{a}}/r_{\textrm{a}}$ is called a shape factor, or sometimes a 'field factor'. The form of equation (13) involving k_a is more commonly used in field ion emission contexts, but is also useful for the needle-like emitters used in electron microscopes and traditional field electron microscopes.

Equation (13) is also sometimes seen in the form

$$\begin{equation} E_{\textrm{a}} = E_{\textrm{free}} / k_{\textrm{a}} , \end{equation} \tag{ 14 }$$

where E_free is the field at the surface of a 'free' sphere charged to the same potential as the needle-like emitter.

For electronically ideal FE systems, empirical values of a characteristic LVCL (usually interpreted as an apex LVCL) can be extracted from current–voltage characteristics as analysed with Fowler–Nordheim plots (e.g. [29]) or (better) with MG plots [23]. However, applying these procedures to non-ideal systems can generate spurious LVCL values [27]. (Also see [18].)

This article is primarily about the 'pure' electrostatics of electronically ideal FE systems.

2.5.1. Emitter geometry.

The main emitter geometry considered in this review is a cylindrically symmetric protrusion, with its apex labelled 'a', having total height h, apex radius of curvature r_a, and defined shape. This protrusion stands on an emitter substrate, with the protrusion symmetry axis passing through a point 'b' on the substrate surface.

The ratio $(h/r_{\textrm{a}})$ has an important theoretical role: we call it the apex sharpness ratio and denote it by σ_a. Thus

$$\begin{equation} \sigma_{\textrm{a}} \equiv h / r_{\textrm{a}} . \end{equation} \tag{ 15 }$$

The substrate plus the protrusion forms the emitter electrode. In the diode geometries discussed in this review, the emitter faces a single counter-electrode, called the collector. In many cases of theoretical interest, the emitter substrate will be one of a pair of well-separated parallel planar plates, of large lateral extent compared to their separation d_sep: as already indicated, this will be called PPP geometry. Usually, it will be assumed that the protrusion height $h \ll d_{\textrm{sep}}$.

However, a later section formally extends analysis to the situation where the protrusion stands at the apex 'b' of a substrate of defined non-planar geometry, with apex radius of curvature much larger than r_a. In this case, the substrate electrostatics can sometimes be the electrostatics of the needle-shaped emitters used in electron microscopes and in field electron and field ion microscopes, i.e. so-called needle geometry. It will be useful to consider this substrate electrostatics first, in section 4.

In principle, field emitter electrostatics may need to consider other geometries—for example those involving concentric spheres, cylindrical wires, blade-shaped emitters or parallel plates of limited lateral extent. However, this review focuses mainly on PPP and needle geometries.

2.5.2. Macroscopic fields.

As discussed in more detail below, when a protrusion is present, then so-called macroscopic-field enhancement factors (MFEFs) can be defined. A local MFEF γ_ML applies to some specific location 'L' on the protrusion surface and is defined by an equation of the form

$$\begin{equation} \gamma_{\textrm{ML}} = E_{\textrm{L}} / E_{\textrm{M}}. \end{equation} \tag{ 16 }$$

Here, E_L is the local electrostatic field at location 'L', and the so-called macroscopic field E_M is an electrostatic field (different in different circumstances) that characterises the FE system geometry more generally.

For the electronically ideal FE systems under consideration, the value of E_M depends on the interelectrode voltage V_M and the system geometry, and is given formally by

$$\begin{equation} E_{\textrm{M}} = -V_{\textrm{M}} / \zeta_{\textrm{M}}, \end{equation} \tag{ 17 }$$

where ζ_M is the (zero-current) macroscopic voltage conversion length. V_M is considered positive if the collector is positive relative to the emitter electrode.

For electronically ideal FE systems, all parameters are defined in the zero-current approximation (i.e. it is assumed there are no significant voltage differences in the system other than that across the HVG). (However, in-system voltage differences can occur in non-ideal systems—see sections 9.1 and 9.2.)

From the above, it follows formally that

$$\begin{equation} \gamma_{\textrm{ML}} = \zeta_{\textrm{M}} / \zeta_{\textrm{L}}.\end{equation} \tag{ 18 }$$

Macroscopic fields fall into two broad types: base-fields (E_B or E_b); and gap-fields (E_G).

A base-field is the field-value that exists at substrate point 'b' as defined earlier, in the absence of the protrusion and of certain effects that would be induced by its presence. (More careful discussions are given later.)

For a given interelectrode voltage, the value of E_B will depend on the functional form and value of the corresponding VCL ζ_B, which in turn will depend on the system geometry. For PPP geometry, ζ_B becomes given by the plate separation d_sep, and E_B becomes the inter-plate field E_P (also called the 'plate field') given by

$$\begin{equation} E_{\textrm{P}} = -V_{\textrm{M}} / d_{\textrm{sep}} \;\; \mathrm{or} \; = -V_{\textrm{P}} / d_{\textrm{sep}} . \end{equation} \tag{ 19 }$$

For notational consistency, in PPP geometry the symbol 'V_M' can be replaced by the symbol 'V_P', although the meanings are the same.

Other geometries involving base-fields are discussed below.

Note a change in terminology as compared with some earlier papers. What was previously called a 'macroscopic field' is now being called an '(inter-)plate field' and is being denoted by E_P rather than E_M (or $|E_{\textrm{P}}|$ would be used instead of F_M), and the term 'macroscopic field' is being used in a more general sense.

A gap-field E_G is defined as the mean field between the apex of the protrusion and an adjacent counter-electrode, between which a voltage V_G exists. E_G is given by

$$\begin{equation} E_{\textrm{G}} = -V_{\textrm{G}} / d_{\textrm{gap}} , \end{equation} \tag{ 20 }$$

where d_gap is the gap length. For electronically ideal systems V_G is equal to the inter-electrode voltage V_M.

Gap-fields are used in various geometrical situations, including close collector adjacency in PPP-type geometry, and more generally in point–point and point–plane geometries.

2.5.3. Macroscopic field enhancement factors.

Associated with the different types of macroscopic field there are different types of MFEF. In some situations different types of MFEF can be defined. For example, in PPP geometry both plate-FEFs and gap-FEFs can be defined. When the protrusion height is small in comparison with both the plate separation d_sep and the gap length d_gap (with the result that d_sep and d_gap have nearly equal values), then the two FEF types have similar values. But when this condition is not true, then the behaviours of the two FEF types diverge (see section 9.8).

More generally, it is often not electrostatically legitimate to compare the numerical values of MFEFs of different kinds, as is sometimes done in the literature.

With field enhancement factors, some considerations and formulae apply to all MFEF types: these common formulae retain the label 'M' in the subscript. Some considerations and formulae apply only to the particular type of MFEF under discussion: in these cases the label 'M' is replaced by an appropriate alternative (e.g. 'P', in the case of base-FEFs in PPP geometry). Some generally applicable formulae are now discussed. Many detailed formulae are given in later sections.

When the characteristic location 'C' is of interest, equations (16) and (18) are replaced by

$$\begin{equation} \gamma_{\textrm{MC}} = E_{\textrm{C}} / E_{\textrm{M}} = \zeta_{\textrm{M}} / \zeta_{\textrm{C}} .\end{equation} \tag{ 21 }$$

When the characteristic location is taken at the emitter apex 'a', then this is replaced by

$$\begin{equation} \gamma_{\textrm{Ma}} = E_{\textrm{a}} / E_{\textrm{M}} = \zeta_{\textrm{M}} / \zeta_{\textrm{a}} .\end{equation} \tag{ 22 }$$

For electronically ideal FE systems, these last two formulae are used for deriving FEF-values from the VCL-values extracted from data-analysis plots, but will yield correct results only if the parameter ζ_M has been determined correctly for the actual geometry of the FE system used.

In general, the value of ζ_M may depend on the whole geometry of the FE system and thus may be difficult to calculate accurately. Most basic modelling and data-analysis use some sort of simple 'special-case' geometry.

In the special case of PPP geometry, the formula for the 'apex plate-FEF (PFEF)' is

$$\begin{equation} \gamma_{\textrm{Pa}} = E_{\textrm{a}} / E_{\textrm{P}} = d_{\textrm{sep}} / \zeta_{\textrm{a}} .\end{equation} \tag{ 23 }$$

In PPP geometry (and also in some other geometries), the formula for the apex gap-FEF (GFEF) is

$$\begin{equation} \gamma_{\textrm{Ga}} = E_{\textrm{a}} / E_{\textrm{G}} = d_{\textrm{gap}} / \zeta_{\textrm{a}} .\end{equation} \tag{ 24 }$$

A merit of PPP geometry is that, provided d_sep is much greater than the protrusion height h (say, $d_{\textrm{sep}}/h \gt 5$) then the values of γ_PL (and hence γ_Pa) are effectively independent of d_sep, and depend only on the protrusion shape. The value of γ_Pa thus characterises the sharpness of the emitter. This basic theoretical situation, and related situations where arrays of identical emitters on the emitter plate are involved, has attracted much theoretical attention, and is a primary focus of this review.

A related situation of practical interest (e.g. [30]) is the concentric cylinders case, where (like PPP geometry) an exact analytical expression can be given for ζ_M.

Usually, it will be difficult to relate a gap FEF to the true electrostatics of the real geometry, and (except for very small gaps) the value of the gap FEF will vary with the gap length (see section 9.8). Thus, gap FEFs do not depend only on the emitter shape, and can have only limited value as emitter characterization parameters.

We stress that formulae (21)–(24) are 'purely electrostatic' formulae that have been formulated for electronically ideal FE systems, and should not be applied blindly to systems that have not been subjected to some form of validity check that the related FE system is ideal.

Figure 1 provides a visual reminder of the hierarchy of terms currently being used to label fields and field enhancement factors.

Zoom In Zoom Out Reset image size

The classical theory of electricity is built on work first comprehensively summarised by Maxwell (see [31]), and on many related later classical developments. However, in the early decades of the 20th century, it became clear that the theory of the behaviour of charged metal and metal-like conductors is (at a more basic level) the theory of electron behaviour in these conductors, and is thus governed by the principles of chemical thermodynamics, quantum mechanics and quantum statistical mechanics.

In particular, as shown by Fowler [32], an appropriately defined electron Fermi level acts as the chemical potential for electrons. In a connected system of electrical conductors in which there are no current-driving elements such as batteries, the condition for static electrical equilibrium is that the Fermi level be constant throughout the system. This is the situation in which the relevant Gibbs function (that includes electron energy and work terms) has its minimum value. Constancy of the Fermi level throughout the system is the atomic-level equivalent of the classical condition that voltage be constant throughout the system.

It has long been known empirically that, when a system is in static electrical equilibrium, the electrostatic field is highest over the most sharply pointed features of the system, with the field being particularly high at the apex of any spike-like feature of the system surface. The fundamental underlying reasons for this are poorly explained both in the electrostatics literature as a whole and in most electricity and magnetism textbooks. At the classical level, the most careful discussion we know of is the 1886 discussion of Robin [33], who (in effect) constructed and minimised an electrostatic energy functional involving the charge distribution across the surface of an irregularly shaped body. A modified version of this approach (using point charges at surface atom positions) has been developed by Rolland et al [34] in the context of atom probe microscopy, in order to model the process of FEV.

In an atomic-level approach, the requirement (as just indicated) is that the relevant Gibbs function have its minimum value and hence that the Fermi level be constant throughout the emitter electrode. When the simplifying assumption is made that the local work function has the same uniform value for all conducting surfaces in the system, then the condition for static electrical equilibrium of the emitter electrode becomes that the electrostatic potential is uniform across the surface of the emitter electrode. Across the surface of the counter-electrode a different constant electrostatic potential can be assumed. These conditions allow the field distribution across the emitter electrode (or relevant parts of it) to be determined using classical electrostatics.

Although there are a few configurations that can be exactly solved analytically (and hence have an important role in validating approximate and numerical solutions), most geometrical configurations cannot be exactly solved analytically. At a classical level, there are two main ways of proceeding: (a) charge-distribution methods; and (b) numerical solution of Laplace's equation.

In a charge-distribution method, one determines an arrangement (positions and magnitudes) of point charges and/or dipoles that yields the correct electrostatic-potential value(s) at a relevant number of chosen points on the active conductor surface (in practice, on the protrusion surface). Then, using these charges and dipoles, and Coulomb-law arguments, one can determine the electrostatic field(s) at a relevant point or relevant points on the protrusion surface. In practice, it is often only the electrostatic field at a protrusion apex that is of significant interest. The simpler forms of this method can be easy to apply and can lead to simple formulae that are physically instructive, although not accurate. The more sophisticated forms lead (in the cylindrically symmetric geometries of interest here) to so-called 'line-charge models'.

Several different techniques exist for solving Laplace's equation numerically. These include finite element methods, finite difference methods and boundary element methods. For numerical work, this paper uses finite element methods exclusively. Relevant basic principles and issues are outlined in section 6.

• (a)  
  For coaxial cylinders an exact formula is

  $$\begin{equation} \zeta_{\textrm{b}} = R_{\textrm{cyl}} \ln(1 + d / R_{\textrm{cyl}}), \end{equation} \tag{ 25 }$$

  where R_cyl is the cylinder radius.
• (b)  
  For concentric spheres an exact formula is

  $$\begin{equation} \zeta_{\textrm{b}} = R_{\textrm{b}} d / (R_{\textrm{b}} + d). \end{equation} \tag{ 26 }$$

  This formula has two limiting cases, as follows.
• (c)  
  When $R_{\textrm{b}} \gg d$ then

  $$\begin{equation} \zeta_{\textrm{b}} \approx d. \end{equation} \tag{ 27 }$$

  In the limit of very large R_b this becomes PPP geometry.
• (d)  
  When $d \gg R_{\textrm{b}}$ then

  $$\begin{equation} \zeta_{\textrm{b}} \approx R_{\textrm{b}} .\end{equation} \tag{ 28 }$$

  In the limit of very large d this becomes free-sphere geometry.As already noted, the practical situation where—as in a traditional field electron or field ion microscope—a needle-shaped emitter faces a distant counter-electrode is called here 'needle geometry'. With needle geometry, it is more convenient to use the form $\zeta_{\textrm{b}} = k_{\textrm{b}}R_{\textrm{b}}$, with $k_{\textrm{b}}\gt1$ and to specify a value for $k_{\textrm{b}}$. The electrostatic argument is that, for a given contribution to apex potential, the contribution to apex field is less when a charge element is distant from the apex. Expressions and values for $k_{\textrm{b}}$ depend on the needle-model used. The commonly used models are as follows.
• (e)  
  Hyperboloid model. This model appears to have been introduced by Eyring et al [35] in 1928 but the Coelho and Debeau 1971 discussion [36] is mathematically clearer. For a hyperboloid facing a distant plane, a distance d from its apex we have

  $$\begin{equation} {\zeta_{\textrm{b}}(\textrm{hyperb})} = k_{\textrm{b}} R_{\textrm{b}} \approx R_{\textrm{b}} \times \tfrac{1}{2} \ln(4d / R_{\textrm{b}}).\end{equation} \tag{ 29 }$$

  This is a satisfactory approximation if $d \gg R_{\textrm{b}}$ [36].
• (f)  
  Paraboloid model. This model appears to have been introduced by Becker [37] in 1951, but the Rose 1956 discussion (see appendix I in [38]) is mathematically clearer. The focal length z_f is the distance from the paraboloid focus to its vertex (i.e. the paraboloid apex), and the apex radius of curvature $R_{\textrm{b}} = 2 z_{\textrm{f}}$. The distance d between the apexes of two confocal paraboloids with focal lengths z_fe (modelling the emitter) and z_fc (modelling the counter-electrode) thus is $d = z_{\textrm{fc}}-z_{\textrm{fe}}$. From equation (15) in [38] it follows that

  $$\begin{equation} {\zeta_{\textrm{b}} = z_{\textrm{fe}} \ln{(z_{\textrm{fc}} / z_{\textrm{fe}})}} = {R_{\textrm{b}} \times \tfrac{1}{2} \ln{[(d/z_{\textrm{fe}})+1]}} , \end{equation} \tag{ 30 }$$

  $$\begin{equation} {\zeta_{\textrm{b}}(\textrm{parab})} = k_{\textrm{b}} R_{\textrm{b}} \approx R_{\textrm{b}} \times \tfrac{1}{2} \ln(2d / R_{\textrm{b}}). \end{equation} \tag{ 31 }$$

  This paraboloid formula is sometimes given incorrectly in the literature.There is a question of what constitutes a good model for a field emitter needle. Probably, a good model should ideally provide adequate representations of all of the following:

  • 1.  
    The local shape of the needle close to its apex.
  • 2.  
    The shape of the needle distant from its apex.
  • 3.  
    The shape (and location) of the counter-electrode in front of the needle-like apex.
  • 4.  
    The shape and location of the counter-electrode 'around the sides of' the needle, both above and below the needle apex.
• (g)  
  Thus, possibly the best analytical model developed so far for fitting the shapes of real needle-like emitters is the sphere-on-orthogonal-cone (SOC) model, introduced by Dyke et al [39]. The slightly complicated details are well described in the Dyke et al paper and will not be reproduced here. This model provides a mathematical description of how (the reciprocal of) the LVCL varies across the emitter surface and down the shank.

Notwithstanding the above analytical formulae, simpler (approximate) numerical formulae have been developed for use in needle geometry, and have found much empirical use.

• (a)  
  The simplest is the well-known 'rough' formula

  $$\begin{equation} k_{\textrm{b}} \approx 5 . \end{equation} \tag{ 32 }$$

  This is stated by Gomer in a 1955 review article [40] (see p 97), but without detailed justification.Using the hyperboloid formula (29), Coelho and Debeau [36] show that k_b varies from 3, for $d/R_{\textrm{b}} = 10^{\,2}$, to 5 for $d/R_{\textrm{b}} = 10^4$, to 6.5 for $d/R_{\textrm{b}} = 10^5$. For a field electron microscope, one might have $d \approx 10^{-2}$ m, $R_{\textrm{b}} \approx 10^{-6}$ m, i.e. $d/R_{\textrm{b}} = 10^4$; thus, $k_{\textrm{b}} \approx 5$ looks a plausible value.
• (b)  
  However, other researchers have suggested that their emitters conform better to shape-factor values that lie somewhere in the range

  $$\begin{equation} 3 \lt k_{\textrm{b}} \lt 8 . \end{equation} \tag{ 33 }$$
• (c)  
  An empirical formula found by Charbonnier ([41], see p 45), apparently by fitting SOC-model profiles to electron micrographs of field electron emitters, is

  $$\begin{equation} k_{\textrm{b}} \approx 0.59 \alpha_{\textrm{c}}^{1/3} (d/R_{\textrm{b}})^{0.13} , \end{equation} \tag{ 34 }$$

  where α_c is the half-angle of the cone that defines the shape of the needle shank. This approach is said to generate shape-factor values close to 5 in many practical cases.
• (d)  
  Biswas [42, 43] has presented a line-charge analysis of needle geometries similar to those discussed above.

We emphasize that the simple formulae discussed above do not apply well to the post-on-plane shapes discussed next.

A particularly simple physical situation is the hemisphere of radius r_H on an infinite plane, in a uniform applied field E_P. This is modelled by placing, at the centre of the hemisphere, a dipole of ISQ polarisability $\alpha_{\textrm{dip}} = 4\pi \epsilon_0 r_{\textrm{H}}^3$. At the hemisphere apex, the field due to the dipole is $E_{\textrm{dip}} = 2E_{\textrm{P}}$. This adds to the applied field, and results in the apex-FEF value $\gamma_{\textrm{Pa}} = 3$. Fields and potentials at other points on the hemisphere can be obtained by combining results due to the usual dipole formulae with those related to the applied field. Note that the total field at any point where the hemisphere intersects the plane is zero, and that at all points on its surface the electrostatic potential of the hemisphere is (necessarily) equal to that of the plane.

Another analytically solvable situation is the cylindrically symmetric hemi-ellipsoid on an infinite plane, in a uniform applied field. Let the lengths of the semi-major axis (or 'height') be h, and of the semi-minor axis (or 'base-radius') be ρ, with $h\gt\rho$. An expression for γ_Pa, in the form used by Forbes et al [46], is as follows. The axis ratio (or aspect ratio) ν of the ellipsoid, a parameter ξ, the apex radius of curvature r_a, and the apex plate-FEF γ_Pa are given by

$$\begin{equation} \nu \equiv \frac{h}{\rho}, \end{equation} \tag{ 35 }$$

$$\begin{equation} \xi \equiv \left(\nu^{\,2} -1\right)^{1/2}, \end{equation} \tag{ 36 }$$

$$\begin{equation} r_{\textrm{a}} = \frac{\rho^{\,2}}{h} = \frac{\rho}{\nu}, \end{equation} \tag{ 37 }$$

$$\begin{equation} \gamma_{\textrm{Pa}} = \frac{\xi^3}{\left\{\nu\ln{\left[\nu + \xi\right]} - \xi\right\} }. \end{equation} \tag{ 38 }$$

It also follows from (35) and (37) that

$$\begin{equation} \nu \equiv \frac{h}{\rho} = \left(\frac{h}{r_{\textrm{a}}}\right)^{1/2} = \sigma_{\textrm{a}}^{1/2}. \end{equation} \tag{ 39 }$$

These results are based on the form given by Latham [47] and attributed to Rohrbach [48]. A more complicated (but equivalent) form of the same result is proved by Kosmahl, equation (10) in [49], starting from the electrostatics of ellipsoids in uniform fields, as given in standard textbooks.

Obviously, in the limit that $\nu \gg 1$ (i.e. for a long thin ellipsoid), then $\zeta \approx \nu$ and equation (38) reduces to the simpler form

$$\begin{equation} \gamma_{\textrm{Pa}} \approx \frac{\nu^{\,2}}{\left[\ln{\left(2\nu \right)} - 1\right] }.\end{equation} \tag{ 40 }$$

Further, since $\nu^{\,2} = h/r_{\textrm{a}}$, result (40) can also be put in the instructive form

$$\begin{equation} \gamma_{\textrm{Pa}} \approx c_{\textrm{HEP}} \times ( h / r_{\textrm{a}} ), \end{equation} \tag{ 41 }$$

where the correction factor c_HEP is given by

$$\begin{equation} c_{\textrm{HEP}} = \frac{1}{\tfrac{1}{2} \ln{\left[ 4\left( h/r_{\textrm{a}}\right) \right]} - 1}.\end{equation} \tag{ 42 }$$

For example, for $h/r_{\textrm{a}} = 100$, $c_{\textrm{HEP}} \approx 0.501$. An alternative approach, based on the use of a non-uniform axial line-charge, has been developed by Pogorelov et al [50, 51]. This may not be an exact general treatment, but leads again to formulae (41) and (42) in the case that $\nu \gg 1$.

On moving away from physical models with exact analytical solutions, one may need to distinguish between: (a) exact physical models; (b) approximate physical models; and (c) related electrostatic models. For the HCP physical model discussed below, a useful approximate physical model is the so-called 'floating sphere at emitter-plane potential' (FSEPP) model (e.g. [52]); a related basic electrostatic model uses only a point charge q placed at the sphere centre.

A slightly generalised version of this basic model is more useful, and is presented here. This generalised model does not require the substrate to be a plane and can be called the basic floating sphere at substrate potential (BFSSP) model.

Consider the situation that exists in the absence of the protrusion (or a model for it), and in the absence of effects due to 'vertically related' electrode images or equivalent (i.e. images in the substrate or in the collector, or equivalent effects), and normally in the absence of any other electrostatic effects would be induced by the protrusion if present (such as depolarization of neighbouring protrusions in a cluster or array, or 'lateral' mirror images in the sides of the simulation box).

Let $(\Delta {\varPhi})_{\textrm{XD}}$ [${\equiv} ({\varPhi}_{\textrm{a}} - {\varPhi}_{\textrm{b}})_{\textrm{XD}}$] be the electrostatic potential difference (between the positions of the protrusion apex 'a' and base 'b') that is generated by 'external' electrostatic sources (i.e. sources other than those directly associated with the protrusion model and its 'vertical' electrode images). Further, let E_Xa be the field generated at 'a' by these external electrostatic sources. We refer to $(\Delta {\varPhi})_{\textrm{XD}}$ as the driving potential difference, and to E_Xa as the external field contribution to the apex electrostatic field.

When the protrusion is present, the point charge in the basic electrostatic model generates a potential difference $(\Delta {\varPhi})_{\textrm{model}}$ between 'a' and 'b'. Thus, the condition that the protrusion/post apex (and hence the sphere apex) be at substrate potential requires that

$$\begin{equation} (\Delta {\varPhi})_{\textrm{model}} + (\Delta {\varPhi})_{\textrm{XD}} = 0 .\end{equation} \tag{ 43 }$$

The solution of this equation generates a value for q and hence a primary contribution E_q to the total electrostatic field $(E_{\textrm{a}})_{\textrm{tot}}$ at the model apex. Thus, in this BFSSP model, the total apex field is

$$\begin{equation} (E_{\textrm{a}})_{\textrm{tot}} = E_q + E_{\textrm{Xa}} , \end{equation} \tag{ 44 }$$

and the apex value $(\gamma_{\textrm{Ba}})_{\textrm{model}}$ of the model-predicted base-field enhancement factor (BFEF) is given by

$$\begin{equation} (\gamma_{\textrm{Ba}})_{\textrm{model}} = (E_{\textrm{a}})_{\textrm{tot}} / E_{\textrm{B}} , \end{equation} \tag{ 45 }$$

where E_B is the base-field, defined here (as before) as given by the field at point 'b', in the absence of the protrusion and of vertically-related electrode images (or equivalent).

In the definition of $(\Delta {\varPhi})_{\textrm{model}}$, the related potentials ${\varPhi}_{\textrm{a}}$ and ${\varPhi}_{\textrm{b}}$ can be taken as measured relative to infinity. In circumstances where $|{\varPhi}_{\textrm{b}}| \ll |{\varPhi}_{\textrm{a}}|$ (which requires $(h/r_{\textrm{a}}) \gg 1$ (say, greater than 10), the above theory can be developed further. In this case it follows from Coulomb's Law and equation (43) that $E_q \approx - (\Delta {\varPhi})_{\textrm{XD}} / r_{\textrm{a}}$, and hence from equation (44) that

$$\begin{equation} (E_{\textrm{a}})_{\textrm{tot}} \approx - (\Delta {\varPhi})_{\textrm{XD}} / r_{\textrm{a}} + E_{\textrm{Xa}} .\end{equation} \tag{ 46 }$$

In many practical circumstances (a possible exception being when emitters in a cluster are close together), the second term on the r.h.s. of equation (46) is much smaller than the first term, and can be neglected. The formula for the model BFEF then simplifies to

$$\begin{equation} (\gamma_{\textrm{Ba}})_{\textrm{model}} \approx - (\Delta {\varPhi})_{\textrm{XD}} / E_{\textrm{B}} r_{\textrm{a}} .\end{equation} \tag{ 47 }$$

As shown below, this relation between the apex BFEF and the driving PD is a key formula that can be applied in various contexts. It is called here the driving potential-difference formula, and applies to both planar and non-planar substrates.

Formula (47) is not an accurate formula, but is considered to display the physics clearly and to provide a 'reasonable first estimate'.

The well-known basic FSEPP-model formula is obtained from equation (47) by replacing E_B by E_P and setting $(\Delta {\varPhi})_{\textrm{XD}} = - E_{\textrm{P}} h$, which yields

$$\begin{equation} {\gamma_{\textrm{Pa}}(\textrm{FSEPP}}) \approx h/r_{\textrm{a}}.\end{equation} \tag{ 48 }$$

The FSEPP model was first applied to field emission by Vibrans [53, 54]. Subsequently, many researchers [46, 48, 52, 55–57] have explored both the basic model and more sophisticated versions that involve placing a dipole at the sphere centre (which makes the sphere surface an approximate equipotential), and/or various image configurations. Review-type discussions are presented in [46, 52, 57].

In earlier stages of the development of FE electrostatics, there was a tendency to treat the FSEPP-derived formulae as useful (albeit approximate) working formulae, and attention was given to getting the small correction terms right. With the development of the more accurate approaches discussed later, this role has fallen away, but formulae derived from equation (47) remain highly useful as models that give insight into the basic physics of sharply pointed conductors—which is not widely understood in the general scientific community.

The two-connected-spheres (TCS) model is another more general version of the FSEPP model and allows basic physical discussion of the electrostatics of protrusions on spherical and on needle-geometry substrates. The substrate is treated as a sphere of large radius R_b that can be modelled electrostatically by placing a large point-charge Q at its centre. From the electrostatics of spheres, it is readily shown that the driving PD now becomes

$$\begin{align} {(\Delta {\varPhi})_{\textrm{XD}}(\textrm{TCS})} & = k_{\textrm{C}} Q [1/(R_{\textrm{b}}+h) - 1/R_{\textrm{b}}] \end{align} \tag{ 49 }$$

$$\begin{align} &\approx (k_{\textrm{C}} Q / R_{\textrm{b}}) \left[ (1+h/R_{\textrm{b}})^{-1} - 1\right] \end{align} \tag{ 50 }$$

where $k_{\textrm{C}} \equiv 1/4\pi \epsilon_0$. On using binomial expansion, and then noting that (in this mode) $(k_{\textrm{C}} Q / R_{\textrm{b}}^{\,2}) = E_{\textrm{B}}$ the above equations yield

$$\begin{align} {(\Delta {\varPhi})_{\textrm{XD}}(\textrm{TCS})} & \approx - E_{\textrm{B}} h , \end{align} \tag{ 51 }$$

$$\begin{align} {\gamma_{\textrm{Ba}}(\textrm{TCS}}) & \approx h/r_{\textrm{a}}. \end{align} \tag{ 52 }$$

That is, the basic TCS result is the same as the basic FSEPP result. (However, as will be shown later, results are different if second-order terms are included.)

In the context of FE from a mercury whisker grown on what appears to be a nearly flat tungsten substrate [58, 59], a version of the TCS model was used by Gomer [58] in 1957 to derive a formula for the shape factor k_a that appears in equation (14) above. In our notation, Gomer's equation (1) becomes

$$\begin{equation} E_{\textrm{a}}/E_{\textrm{free}} = 1/k_{\textrm{a}} = h/(R_{\textrm{b}}+h) \approx h/R_{\textrm{b}}. \end{equation} \tag{ 53 }$$

Since, for two spheres at the same potential, $E_{\textrm{free}} / E_{\textrm{b}} = R_{\textrm{b}} / r_{\textrm{a}}$, equation (53) is equivalent to equation (52).

A different approach to the mathematical analysis of the two-connected-sphere situation has been presented by Feynman et al [60] and also by Lewin [61]. However, their formula predicts the apex FEF for the smaller sphere to become infinite as R_b becomes large; thus, we regard their analyses as either physically inappropriate or physically incorrect [62]. This point has previously been made by Fricker [63].

A feature learnt by comparing the Feynman treatment with the treatment here is the following. For a given charge under analysis, best practice is to write down expressions for the potential difference between two well-defined points in the system model (even if this results in a small term that is later discarded); it can be poor practice to write expressions for the 'potential relative to infinity'. This is because the latter practice can lead to the accidental omission of apparently small terms that turn out to be important in some contexts.

It is clear that the apex sharpness ratio $h/r_{\textrm{a}}$ is the most important factor in the electrostatics of sharp points and posts. Since equations (48) and (52) are not accurate, it is convenient to introduce a correction factor c_a and write

$$\begin{equation} {\gamma_{\textrm{Pa}}(\textrm{CP})} = c_{\textrm{a}} (h/r_{\textrm{a}}). \end{equation} \tag{ 54 }$$

This has been called the conducting post formula.

When a conducting but electrically isolated 'floating cylindrical rod' of length L_rod is situated in an uniform applied field E_P and aligned parallel to it, then a planar equipotential surface passes through the centre of the rod, normal to its length. The top half of this configuration is equivalent to conducting post geometry. Hence, the apex FEF $\gamma_{\textrm{Pa}}(\textrm{FR})$ for the floating rod is given by

$$\begin{equation} {\gamma_{\textrm{Pa}}(\textrm{FR})} = \tfrac{1}{2} c_{\textrm{a}} (L_{\textrm{rod}}/r_{\textrm{a}}). \end{equation} \tag{ 55 }$$

This has been called floating rod formula.

For reasons that are not entirely clear, these formulae do not normally appear in electricity and magnetism textbooks, and appear not to be widely known. They promise to have significant applications outside FE when discussing electrical phenomena (particularly breakdown phenomena) associated with pointed conductors, even in the approximation $c_{\textrm{a}} = 1$.

From about 2000, serious attempts were made to develop more accurate analyses of the HCP model. These included a complicated charge-distribution model devised by Kokkorakis et al [64], and FEM treatments by Edgcombe and Valdrè (EV) [65–67]. For a review in 2003 and a comparison of detailed formulae, see [46]. An useful conclusion [46, 65] was that the approximation

$$\begin{equation} c_{\textrm{a}} \approx 0.7\end{equation} \tag{ 56 }$$

was valid to within about 25% in the parameter range $30 \lt (h/r_{\textrm{a}}) \lt 2000$.

In PPP geometry, a related formula for the apex shape factor k_a can be obtained from the conducting post formula, by making use of equation (23) in the form

$$\begin{equation} k_{\textrm{a}} = \frac{d_{\textrm{sep}}}{\gamma_{\textrm{Ma}} r_{\textrm{a}} }. \end{equation} \tag{ 57 }$$

The result is

$$\begin{equation} {k_{\textrm{a}}(\textrm{PPP})} \approx \frac{d_{\textrm{sep}}}{c_{\textrm{a}} h }. \end{equation} \tag{ 58 }$$

As illustrative values, we can take $d_{\textrm{sep}} = 25\,\mu$m, $h = 1\,\mu$m, $c_{\textrm{a}} = 0.7$, which yields $k_{\textrm{a}}(\textrm{PPP}) \approx 36$. This shows that the needle-substrate approximation $k_{\textrm{b}} \approx 5$ does not work for a post in PPP geometry. The failure of this approximation to describe a complex tip geometry was previously shown by EV [68].

For the hemisphere-on-plane (HSP) model, figure 3 illustrates how the calculated apex-PFEF value γ_Pa depends on the simulation domain dimensions. FEM simulations need to use local densities of mesh elements sufficient to ensure that related numerical errors are smaller than errors due to the finite domain dimensions. Greater densities of FEM elements need to be used where the electrostatic fields are larger. Near the HSP emitter apex the use of element sizes smaller than $r_{\textrm{a}}/30$ is usually enough to evaluate electrostatic fields with the desired precision (see inset of the figure 3); this condition (but written as 'size $\lt r_{\textrm{a}}/30$') also applies to the apexes of other emitter shapes.

Zoom In Zoom Out Reset image size

The total percentage error (T%E) ε_t in the numerical apex-PFEF value $(\gamma_{\textrm{Pa}}^{\textrm{num}}$), relative to the analytical apex-PFEF value ($\gamma_{\textrm{Pa}}^{\textrm{anal}}$), can be defined by the formula (which ensures that T%E is usually positive):

$$\begin{equation} \epsilon_{\textrm{t}} \equiv \frac{(\gamma_{\textrm{Pa}}^{\textrm{anal}} - \gamma_{\textrm{Pa}}^{\textrm{num}}) } {\gamma_{\textrm{Pa}}^{\textrm{anal}}} \times 100\%. \end{equation} \tag{ 59 }$$

In this review, errors in numerical FEF values will usually be defined in this general way.

For the HSP model the correct apex-FEF is exactly 3 (see section 5.1). It is found that both the radius A and the height B of the simulation domain influence the total percentage error in the simulated apex-FEF value, with the T%E decreasing significantly as the box dimensions increase. If one of the box dimensions (say the height B) is taken as large, and thus does not contribute significantly to the T%E, then it is found that for moderate values of the other box dimension (in this case the radius A), the contribution ε_A of this dimension (as defined by the value of ε_t in these circumstances) to the T%E falls off with a power law close to 3. The same effect is found if the roles of box height and radius are reversed. These results are illustrated in figure 4, where the plots are made against the (dimensionless) normalised quantities $A/r_{\textrm{H}}$ and $B/r_{\textrm{H}}$, where r_H is the hemisphere radius.

Zoom In Zoom Out Reset image size

For the lateral boundary, the origin of the power-law exponent 3 can be understood by considering a small angular element of the circular mirror-dipole. This acts as in the same way as the distant dipole in the two-post-like-emitter situations analysed in [72, 73]. The apex-field reduction for the chosen emitter shape is related to the value of the depolarizing field at point 'b' due to the mirror-dipole; for sufficiently large values of box-radius A, this depolarising dipole-field falls off as $(A^{-3})$. Each element of the circular mirror-dipole behaves in this way, so the total depolarizing field behaves in this way. This argument is very general, so all post-like shapes are expected to have $(A^{-3})$-type limiting behaviour.

For the top boundary, the error behaviour is also an effect associated with a dipole inverse-cube law, but is of a different kind. At sufficient large distances in any direction, the additional field due to a post-like emitter (above that due to the uniform applied field in which it is situated, in PPP geometry) can be modelled as due to a point dipole at 'b'. (With a hemisphere on a plane, the modelling is exact.) Hence, at the top boundary, the error in field specification created by imposing a uniform field at the boundary equal to the applied macroscopic field is approximately given by the field distribution on the boundary due to a dipole at 'b'.

It appears that in the simulations this error translates into an error in the apex field value that falls off with the inverse cube of box height. The exact electrostatic mechanism by which this simulation error is created when solving Laplace's equation is not yet clear, but the observed inverse-cube dependence on box height looks entirely plausible physically. Again, the argument is very general, so all post-like shapes are expected to have $(B^{-3})$-type limiting behaviour.

The data in figure 4 can usefully be presented in another way, which better illustrates the combined effect of the two sources of error. This is done in figure 5. For several values of $A/r_{\textrm{H}}$, the T%E ε_t is plotted as a function of $B/r_{\textrm{H}}$. The left-hand (sloping) part of these curves shows the behaviour of the partial percentage error ε_B , as already plotted in figure 4. However, when this error ε_B falls below the partial percentage error ε_A for the chosen value of $A/r_{\textrm{H}}$, the T%E ε_t 'stabilises' at the relevant value of ε_A . The critical value $(B/r_{\textrm{H}})_{\textrm{c}}$ at which this occurs (defined as the intersection of the extrapolated horizontal line with the slope) is shown in figure 6 as a function of $A/r_{\textrm{H}}$, and is well modelled by the relation $B/r_{\textrm{H}} = 0.67 A/r_{\textrm{H}}$. This suggests that this critical value corresponds to a critical domain height/radius ratio of $(B/A)_{\textrm{c}} = 0.67$, valid for all values of A and B likely to be used in practice. Equivalent plots could be drawn, and equivalent arguments could be made, with the roles of domain height and radius interchanged. In the regime where ε_A and ε_B are approximately equal, the two errors combine to give the total error ε_t, but ε_t is not the simple sum of ε_A and ε_B .

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the purposes of efficient simulation, the need is to decide a maximum tolerable T%E, denoted here by ε_tol, and to determine corresponding minimum values of $A/r_{\textrm{H}}$ and $B/r_{\textrm{H}}$. Figure 5 shows that the curve for $A/r_{\textrm{H}} = 11$ corresponds to $\epsilon_{\textrm{t}} = 0.12\%$, for $B/r_{\textrm{H}}$ greater than about 10. However, for ε_t-values that do not correspond to the values of $A/r_{\textrm{H}}$ used in figure 5, a more general approach is needed. This is achieved as follows. It is a reasonable hypothesis that the point at which ε_A and ε_B are equal is determined by a particular value of $B/A$ somewhere around unity; it is also known that for the 'stabilised' part of the curves in figure 5 the percentage error ε_A falls off as $(A/r_{\textrm{H}})^{-3}$. These things suggest using the data from figure 5 to make a trial plot of $\epsilon_{\textrm{t}} \times (A/r_{\textrm{H}})^3$ versus $B/A$. As figure 7 shows, it is found that all the curves in figure 5 then collapse into a single 'slightly broadened curve' applicable to any value of $(A/r_{\textrm{H}})$.

Zoom In Zoom Out Reset image size

The highest of the horizontal sections (or 'plateaus') of these re-plotted curves corresponds to the value $\epsilon_{\textrm{t}} \times (A/r_{\textrm{H}})^3 = 175$. Hence, in a first approach one can use the formula that the minimum value $(A/r_{\textrm{H}})_{\textrm{min}} = (175/\epsilon_{\textrm{tol}})^{1/3} \approx 5.59/(\epsilon_{\textrm{tol}})^{1/3}$. The minimum value of $B/A$ is determined by a point at the left-hand end of the plateau, and figure 6 suggests that a first estimate of this should be taken as $(B/A)_{\textrm{min}} = 0.67$ (or perhaps a little higher). In fact, this value corresponds to the intersection of the sloping line and the extrapolated horizontal lines; in order to be on the truly horizontal part of the broadened curve, a slightly higher value needs to be taken. In practice, we use the value $(B/A)_{\textrm{min}} = 0.83$ derived in [71], by a slightly different argument.

It follows that the minimum value $(B/A)_{\textrm{min}}$ would be given by $(B/A)_{\textrm{min}} \approx 4.63/(\epsilon_{\textrm{tol}})^{1/3}$. It has been checked by numerical simulations that when the values of A and B are chosen according to these criteria, then the actual T%E ε_t is less than the chosen value of ε_tol. However, the view is taken here that using three significant figures in these formulae is 'overkill', so we have rounded them up to the nearest integer, and (for the hemisphere-on-plane model) recommend here the MDD formulae

$$\begin{equation} \left(\frac{A}{r_{\textrm{H}}}\right)_{\textrm{MDD-HSP}} = \frac{6}{\sqrt[3]{\epsilon_{\textrm{tol}}}}, \end{equation} \tag{ 60 }$$

$$\begin{equation} \left(\frac{B}{r_{\textrm{H}}}\right)_{\textrm{MDD-HSP}} = \frac{5}{\sqrt[3]{\epsilon_\mathrm{tol}}}.\end{equation} \tag{ 61 }$$

The above treatment is based on calculating apex-PFEFs. As illustrated in figure 8, which shows the numerically calculated total percentage error ε_t as a function of polar angle θ, we have demonstrated that these MDD values will also guarantee that PFEF-values at positions down the sides of the emitter meet the tolerance used in the equations. (Except that poorer precision is achieved for angles close to $\pi/2$ rad, where the PFEF-value becomes close to zero.)

Zoom In Zoom Out Reset image size

In a HCP model, the apex radius of curvature is equal to the radius of the cylinder and thus to the radius of the protrusion base. This may be a good model for some emitters, but when the emitter shape is such that the apex radius is significantly less than the base radius, then this model does not work well. In this case, the HEP emitter model discussed in this section may be a better choice.

As illustrated in figure 9(a), the HEP-model geometry is obtained by rotating a quarter of an ellipse around its longer axis, to yield the shape shown in figure 9(b). This shape is characterised by the half-axis lengths of the ellipse, which become the emitter base-radius ρ and total height h. For MDD discussions, this shape is usefully characterised by its aspect ratio $\nu \equiv h/\rho$ (see section 5.2). Equation (38) earlier provides an expression for the corresponding apex PFEF.

Zoom In Zoom Out Reset image size

For the HEP emitter, formulae for MDDs have been found by empirical numerical methods, as follows. Simulation-box dimensions are now described by the revised normalised variables $A/h$ and $B/h$. That is, the hemiellipsoid total height h has replaced the hemisphere radius r_H as the normalising length. Obviously, these new quantities reduce to the old ones in the case that ν = 1. In initial trials, for a set of given values of tolerated error ε_tol, minimum dimensions were derived from equations (60) and (61), with the new normalised variables replacing those in the equations. Figure 10 shows how the total error ε_t decreases as the aspect ratio ${\nu }$ increases. As ${\nu }$ increases from 1, the decrease is initially approximately exponential (see inset), but tends to an asymptotic value ε_lim for $\nu \gtrsim 20$ [71].

Zoom In Zoom Out Reset image size

It was shown in [71] that fitting an equation to these results yields the formula

$$\begin{equation} \frac{\epsilon_{\textrm{t}}(\nu)}{\epsilon_{\textrm{t}}(\nu = 1)} \approx 0.2 + 0.8\exp{\left[-0.345\left( \nu - 1\right)\right]}.\end{equation} \tag{ 62 }$$

This result enables us to give the following empirical formulae for the MDDs for the HEP emitter:

$$\begin{align} \left(\frac{A}{h}\right)_{\textrm{MDD-HEP}} = 6\times \sqrt[3]{\frac{0.2 + 0.8\exp{\left[-0.345\left(\nu - 1\right) \right] }}{\epsilon_{\textrm{tol}}}} \end{align} \tag{ 63 }$$

$$\begin{align} \left(\frac{B}{h}\right)_{\textrm{MDD-HEP}} = 5\times \sqrt[3]{\frac{0.2 + 0.8 \exp{\left[-0.345\left(\nu - 1\right)\right]}}{\epsilon_{\textrm{tol}}}}.\end{align} \tag{ 64 }$$

In practice, it is often useful (as in the following section) to have the MDD values expressed as a function of the apex sharpness ratio $\sigma_{\textrm{a}} = h/r_{\textrm{a}}$. From earlier, we have (for the HEP model) $\sigma_{\textrm{a}} = \nu^{\,2}$, so we can define a function $f(\sigma_{\textrm{a}})$ by

$$\begin{equation} f(\sigma_{\textrm{a}}) = 0.2 + 0.8 \exp\left[-0.345\left(\sigma_{\textrm{a}}^{1/2} - 1\right)\right], \end{equation} \tag{ 65 }$$

and the above equations become replaced by

$$\begin{equation} \left(\frac{A}{h}\right)_{\textrm{MDD-HEP}} = 6 \times \sqrt[3]{\frac{f(\sigma_{\textrm{a}})}{\epsilon_{\textrm{tol}}}} \end{equation} \tag{ 66 }$$

$$\begin{equation} \left(\frac{B}{h}\right)_{\textrm{MDD-HEP}} = 5\times \sqrt[3]{\frac{f(\sigma_{\textrm{a}})}{\epsilon_{\textrm{tol}}}}.\end{equation} \tag{ 67 }$$

We now report the results of using these formulae. Figure 11(a) shows a comparison between the exact analytical apex-FEF-value $\gamma_{\textrm{a}}^{\textrm{anal}}$ and numerical values $\gamma_{\textrm{a}}^{\textrm{num}}$ found using equations (65)–(67), for a required ε_tol-value of 0.1%. Figure 11(b) shows the corresponding total percentage errors found in practice. The numerical reasons for the unexpected irregularities in the error graph are not clear, but the T%E is less than 0.1% in all cases in the range of $1.01\lesssim \sigma_{\textrm{a}} \lesssim 1000$. This will be our target precision for other shapes as well.

Zoom In Zoom Out Reset image size

We now consider other single-protrusion shapes. If shapes of this kind, other than the HP and HEP models, are to be simulated, then a method is needed to define MDDs. Our chosen procedure is as follows. As shown in figure 12, a hemiellipsoidal shape is defined that circumscribes the intended emitter shape. The height and apex sharpness ratio of this circumscribing hemiellipsoid are then used in equations (65)–(67) in order to specify MDDs. The thinking behind this approach is that the closer any part of the emitter shape is to the simulation box boundaries, then the worse will be the errors. Hence, since the shape inside the circumscribing hemiellipsoid is everywhere less distant from the lateral simulation-box boundary, the errors are expected to be less than they would be for the hemiellipsoid.

Zoom In Zoom Out Reset image size

It may be that, for some shapes, different procedures could be specified that would lead to slightly smaller MDD values; however, our view is that more detailed investigations of this kind can wait until there is a demonstrable need for them.

As an example, figure 12 shows the much studied HCP emitter, with height h_HCP and apex-radius r_HCP, located within the circumscribing hemiellipsoid. In the case shown, the height h_HEP and apex radius of the hemiellipsoid are taken equal to those of the HCP model emitter. But, in fact, the parameters of any hemiellipsoid that totally encloses the emitter of interest will guarantee the condition $\epsilon_{\textrm{t}} \lt \epsilon_{\textrm{tol}}$.

With emitter systems that are not rotationally symmetric (for example, small emitter clusters) a similar procedure can be applied. One can draw a hemiellipsoid that encompasses the whole emitter system, and use the MDD-HEP formulae to decide working values of $A/h_{\textrm{HEP}}$ and $B/h_{\textrm{HEP}}$ that will implement a given tolerance percentage error. In some cases it may be possible to draw hemiellipsoids of this kind in more than one way, so these working values (although fit for purpose) may not represent true minimum MDD values.

As part of the work described above we found [71], and later used [74], an even more accurate numerical method of finding PFEF values—the so-called extrapolation method. This section illustrates the method and validates it by using the known analytical results for the HEP model, stated earlier.

For a fixed value (50 nm) of apex-radius r_a, and a given value of apex sharpness ratio σ_a, equations (66) and (67) can be used to determine MDDs for a chosen value of the tolerated error ε_tol. Using these domain dimensions then generates a specific prediction for the HEP-model apex-PFEF γ_Pa. The results of carrying out this procedure, for several σ_a-values, and for a large number of ε_tol-values for each σ_a-value, are shown in figure 13. For any chosen σ_a-value, a regression line can be fitted to the related points and the 'extrapolated' PFEF-value that corresponds to $\epsilon_{\textrm{tol}} = 0$ determined.

Zoom In Zoom Out Reset image size

It turns out that this extrapolated numerical value of the apex-PFEF γ_Pa is highly accurate, as is illustrated in table 2. All percentage differences from the known analytical HEP-model results are less than 0.001%, over a wide range of values of apex sharpness ratio σ_a. This finding suggests strongly that this extrapolation methodology should produce reasonably accurate PFEF-values for post-like shapes for which there is no known exact analytical result.

Table 2. Comparisons of apex-PFEF (γ_Pa) values for the hemiellipsoid-on-plane (HEP) model. For specific values of the apex sharpness ratio σ_a, the numerical result from the extrapolation method is compared with the known analytical result given by equations (35)–(39).

σa	γPa—extrapolated	γPa—analytical	% difference
5	6.514 29	6.514 30	$-0.000{\;}15\%$
10	9.816 57	9.816 64	$-0.000{\;}71\%$
25	17.914 36	17.914 41	$-0.000{\;}28\%$
300	117.249 49	117.249 12	$+0.000{\;}32\%$
1000	317.259 44	317.258 71	$+0.000{\;}23\%$
3000	811.1931	811.1916	$+0.000{\;}18\%$

We note that this is a procedure that has been established in PPP geometry. It may be that analogous procedures exist for other geometries, but this is not yet established.

Obviously, the treatments here have employed a cylindrical simulation box. The other main option (which may seem more obvious to some readers) is to use a cuboidal simulation box with a square base. For some problems, in particular the analysis of depolarisation effects in regular square arrays, a cuboidal box of this kind is needed. But for analysis of single protrusions and small clusters the cylindrical box is markedly superior. This is for two main reasons. First, the higher symmetry of the cylindrical box makes meshing problems easier, and very significantly reduces the amount of computing time required to carry out a particular analysis. Second, at least with the finite-element package that we are using, we find in practice that the numerical accuracy and consistency of our results is better when using cylindrical simulation boxes.

As has been shown clearly in this review and in earlier work, the accuracy of electrostatic finite-element analysis can depend significantly on the relative sizes of the object being simulated and the simulation box. For this reason, we recommend that ideally the reported results of electrostatic finite-element simulations should always be accompanied by details of (a) the shape and dimensions of the simulation box, (b) the target accuracy of the simulations, and (c) the methods used to decide on the box dimensions.

Following the pre-2003 work discussed earlier, the HCP model was investigated numerically by Read and Bowring (RB) [75], using a Boundary Element Method, and by several research groups [71, 74, 76–79], using Finite Elements Methods. Resulting fitted formulae for the apex PFEF were presented either as a power series in $(h/r_{\textrm{a}})$ or as an expression of the form

$$\begin{equation} \gamma_{\textrm{Pa}} \approx C \left( D + \frac{h}{r_{\textrm{a}}} \right)^{\kappa}.\end{equation} \tag{ 68 }$$

Edgcombe and Valdrè (EV) [46, 65], using FEM, had reported that, in the range $4\leqslant h/r_{\textrm{a}} \leqslant 3000$, their numerical results for γ_Ma were adequately fitted, to within 3%, by taking C = 1.2, D = 2.15 and κ = 0.90. RB [75], using BEM, claimed that their numerical results would be fitted, with accuracy better than 0.4%, by taking C = 1.0782, D = 4.7 and κ = 0.9152, at least for $31\leqslant h/r_{\textrm{a}} \leqslant 3001$.

Using our MDD extrapolation method, as described in section 7.4, it has been possible to find γ_Pa for the HCP model with high accuracy. A least-squares fitting procedure enables a four-parameter fitting formula, which—following [74]—we prefer to write in the form

$$\begin{equation} \gamma_{\textrm{Pa}} = a\left[b + \left( \frac{h}{r_a} \right)^{c}\right]^d, \end{equation} \tag{ 69 }$$

where a = 0.86391, b = 0.97756, c = 0.52989 and d = 1.77667. In the range $4\leqslant \sigma_{\textrm{a}} \leqslant 3000$, this formula reproduces the extrapolation method result with accuracy better than 0.41%.

Table 3 compares these various numerical and formula-based apex-PFEF estimates, and also provides a value for the correction factor c_a in equation (54). Our conclusion is that our values and the RB values are better than the EV values, and that in nearly all cases formula (69) is slightly better than formula (68). The general consistency of our results and the RB results gives additional confidence that these results for the HCP model are numerically accurate. Figure 14 illustrates how the correction factor c_a depends on the apex sharpness ratio σ_a.

Zoom In Zoom Out Reset image size

The good fitting properties of equation (69) suggests that it could perhaps be applied to other emitter shapes. As an example, we consider a paraboloidal emitter. Figure 15(a) shows numerical and fitted apex-PFEF values, as functions of apex sharpness ratio σ_a. Least-squares fitting using equation (69) yields a = 0.571, b = 1.125, c = 0.42 and d = 2.133. Figure 15(b) plots the total percentage error ε_t, showing that the four-parameter form (69) is a very good fitting function in this case too.

Zoom In Zoom Out Reset image size

A second example is the spherically rounded (truncated) cone (SRC) emitter model, originally discussed in [74]. In this case, the vertex half-angle θ is an additional independent variable, and the parameters in our four-parameter fitting formula become functions of θ, as shown in the inset to figure 16. For 'a' we assume a dependence of the form

$$\begin{equation} a = \omega_{\textrm{a}} + \lambda_{\textrm{a}} \theta^{\mu_{\textrm{a}}}, \end{equation} \tag{ 70 }$$

where ω_a, λ_a and µ_a are fitting parameters. Similar expressions are obtained for 'b', 'c' and 'd'. For fitting over the range $2^{\circ} \leqslant \theta \leqslant 20^{\circ}$, values of the resulting fitting parameters are shown in table 4.

Zoom In Zoom Out Reset image size

Other significant cylindrical-symmetry emitter shapes are considered in the literature but are not discussed here. In particular, these include open-cap nanotubes [80–82] and cusp-like emitter shapes. We expect our MDD arguments to apply to these emitters. They have already been used to discuss open CNTs [82].

When there is significant electrical resistance in the current path between the HVG and the emitter apex, this can cause a Fowler–Nordheim (or related) data plot to bend downwards at the left-hand (high voltage, etc) side. This behaviour is sometimes called 'saturation', and can affect the plot interpretation.

Series resistance effects of this general kind have been discussed theoretically in many papers, including [29, 83–89].

The effect occurs because an electronically ideal FE system has an intrinsic emission resistance R_e that is associated with the emission process itself and is voltage and current dependent.

To illustrate the effect theoretically, we use extended MG FE theory and write an expression for the 'emission current' I_e in terms of the 'emission voltage' V_e (i.e. the voltage applied between the emitting surface and a counter-electrode) in the following simplified way

$$\begin{equation} I_{\textrm{e}} = A_{\textrm{fC}}^{\,\textrm{SN}} C V_{\textrm{e}}^{\,2}\exp[-B/V_{\textrm{e}}], \end{equation} \tag{ 71 }$$

where $A_{\textrm{fC}}^{\,\textrm{SN}}$, B and C will normally be weakly voltage-dependent quantities. The emission resistance R_e is then given by

$$\begin{equation} R_{\textrm{e}} = V_{\textrm{e}}/I_{\textrm{e}} = (V_{\textrm{e}} A_{\textrm{fC}}^{\,\textrm{SN}} C)^{-1} \exp[B/V_{\textrm{e}}] .\end{equation} \tag{ 72 }$$

At lower emission voltages the dominating term in this equation will be the exponent. When V_e is very low, then the exponent and the emission resistance will be very large. As the emission voltage increases, the emission resistance steadily decreases. At some point the decreasing emission resistance may become comparable with the series resistance, and the FN plot ceases to be linear.

The illustrative discussion below is based on PPP geometry, but analogous effects occur in other geometries.

Resistance in series with the emission resistance can arise in several ways. Resistances can be associated with: (a) the post itself; and/or (b) a contact resistance at the point where the post joins the substrate. These are internal series resistances (i.e. between the plates) and can affect the between-plates field distribution, as discussed in section 9.2.

Resistances can also be associated with (for example): (c) a resistive (non-metallic) substrate on which the post-like emitter stands; (d) a series transistor fabricated for current-limiting purposes; (e) poor electrical contacts (and/or the use of resistive materials) in the mechanical emitter support system; and/or (f) a high-value 'safety resistor' (if any voltage drop across this is not subtracted off before processing experimental results). All these constitute external series resistances.

For illustration, we now consider the case of an emission resistance in series with a constant external series resistance R_sr. This will illustrate the principles, but the analysis of real systems could be vastly more complicated.

Even for the simple system just indicated, the algebra is slightly messy, so is not presented here. The results of illustrative simulations are shown in figure 17, as a FN plot using the predicted values of measured current and voltage. It is clear that the plot slope is constant in the low-voltage region, but that its magnitude then decreases as the measured voltage (and hence the emission voltage) increase. At some point the decreasing emission resistance may become comparable with the series resistance, so the FN plot ceases to be linear.

Zoom In Zoom Out Reset image size

For an electronically ideal FE system, or in a low-voltage range of a more complicated system that is acting in an ideal fashion in this low-voltage range, the usual MG/EMG formula for extracting a value of the true characteristic plate-FEF γ_PC is

$$\begin{equation} \gamma_{\textrm{PC}}^{\textrm{extr}}({\textrm{true}}) = - {\textrm{s}_{\textrm{t}}} {b {\phi}^{3/2} d_{\textrm{sep}}} / S_{\textrm{Vm}}^{\textrm{fit}}({\textrm{low}}V_{\textrm{m}}), \end{equation} \tag{ 73 }$$

where $S_{\textrm{Vm}}^{\textrm{fit}}({\textrm{low}}V_{\textrm{m}})$ is the slope of the line fitted to the low-voltage part of the FN plot, and the other symbols have their usual meanings.

Clearly, if the slope of a line fitted to a 'high-voltage' part of the plot were used is this formula, then the formula would yield a false value of the actual PFEF.

In formal principle, one could get a correct formula for interpreting the high-voltage part of the plot by replacing equation (73) by

$$\begin{equation} \gamma_{\textrm{PC}}^{\textrm{extr}}({\textrm{true}}) = - {\sigma_{\textrm{sr,t}}} {b {\phi}^{3/2} d_{\textrm{sep}}} / S_{\textrm{Vm}}^{\textrm{fit}}, \end{equation} \tag{ 74 }$$

where $\sigma_{\textrm{sr,t}}$ is an appropriately defined/determined fitting value of a generalised slope correction factor.

In reality, $\sigma_{\textrm{sr,t}}$ is expected to be a relatively sensitive function of the precise range of V_m-values over which the measurements have been taken, and accurate prediction of the effective correction factor to be used in any particular case looks to be exceptionally difficult.

Also note that experimentally observed 'saturation effects' are not necessarily due to external-series-resistance effects. As discussed in [88], the 'voltage-difference' effects discussed in the next section may also be a possibility, and may in some or many cases provide a more plausible explanation of 'saturation' than external-series-resistance effects.

9.2.1. Background.

9.2.2. Basic PPP-geometry theory.

Discussion in PPP geometry is somewhat simpler. For a resistive post in PPP geometry, the simple formula (79) below was first obtained by Minoux et al [90], by fitting to the results of numerical simulations on a model of a CNT. Later, Forbes [88] presented a simple physical derivation of the formula, and argued that the physical effect was a general phenomenon in finite-current electrostatics. The present treatment is a slightly improved version of the [88] treatment.

When a voltage difference exists between 'a' and 'b', equation (43) is replaced by

$$\begin{equation} (\Delta {\varPhi})_{\textrm{model}} + (\Delta {\varPhi})_{\textrm{XD}} = (\Delta {\varPhi})_{\textrm{d}} = V_{\textrm{d}} , \end{equation} \tag{ 75 }$$

where $(\Delta {\varPhi})_{\textrm{d}}$ [$= {\varPhi}_{\textrm{a}} - {\varPhi}_{\textrm{b}}$] is the electrostatic potential difference (ESPD) corresponding to V_d, with constant work function over the system surfaces being assumed (for simplicity), as before.

For the case considered here of a single emitting protrusion, $(\Delta {\varPhi})_{\textrm{XD}} = - E_{\textrm{P}} h$, as before. Further, if the secondary contribution that E_P makes to the total apex ES field (here denoted simply by E_a) is disregarded, as when deriving equation (44), then—for the basic FSEPP model in use—we have ${(\Delta {\varPhi})_{\textrm{model}} = \Delta {\varPhi}(\textrm{FSEPP)}} {\approx E_{\textrm{a}} r_{\textrm{a}}}$. Equation (75) thus reduces to

$$\begin{equation} E_{\textrm{a}} r_{\textrm{a}} \approx E_{\textrm{P}} h + V_{\textrm{d}} . \end{equation} \tag{ 76 }$$

Algebraic manipulation then yields an expression for the operating apex PFEF:

$$\begin{align} \gamma_{\textrm{Pa}}^{\textrm{op}} \equiv E_{\textrm{a}}/E_{\textrm{P}} \approx (h/r_{\textrm{a}}) (1 + V_{\textrm{d}} / E_{\textrm{P}} h) \end{align} \tag{ 77 }$$

$$\begin{align} \approx \gamma_{\textrm{Pa}}^{\textrm{zc}} \; (1 + V_{\textrm{d}} / E_{\textrm{P}} h), \end{align} \tag{ 78 }$$

where equation (48) has been used to replace $(h/r_{\textrm{a}})$ by the zero-current apex PFEF $\gamma_{\textrm{Pa}}^{\textrm{zc}}$. (This has also been called the 'small-current apex FEF', but it is now thought that 'zero-current' is a better description.)

When this equation is applied to FE, then E_P is negative and V_d is positive. For field ion emitters the signs are reversed. To get a formula that applies to both polarities, equation (77) is re-written as

$$\begin{equation} \gamma_{\textrm{Pa}}^{\textrm{op}} \approx (1 - |V_{\textrm{d}}| / |E_{\textrm{P}}| h) \; \gamma_{\textrm{Pa}}^{\textrm{zc}}.\end{equation} \tag{ 79 }$$

This equation shows clearly that significant apex-PFEF reduction occurs when the voltage change $|V_{\textrm{d}}|$ along the emitter becomes a significant fraction of the magnitude $|E_{\textrm{P}} h|$ of the driving PD (which is provided, in the single-emitter case, by the applied (inter-)plate field alone).

We emphasise that, although the basic theory given here is based on the behaviour of a single post-like protrusion in PPP geometry, this FEF-reduction effect is a general phenomenon in finite-current electrostatics, and should occur with more-or-less any form of resistive emitter with a surface that possesses sharp features.

9.2.3. Comments.

Another effect that can reduce apex-field values (and increase apex-VCL values), for a given voltage difference between the emitter apex and the counter-electrode, is the presence of field-emitted vacuum space charge (FEVSC). In the past, this has been discussed using PPP geometry, with one plate taken as the emitter. Obviously, this is a one-dimensional (1D) model and does not strictly apply to real needle-shaped or post-shaped emitters, but it has merit as an indication of the maximum size of FEVSC effects. In the FEVSC model discussed here, the voltage between the plates is denoted by V_e and termed the emission voltage.

The original work on vacuum space-charge was done in the context of thermally emitted positive ions, and is associated with the names of Child [101] and Langmuir [102, 103]. However, in Child–Langmuir theory the BC at the emitter surface is that the electrostatic field be zero. Thus, strictly, Child–Langmuir theory does not apply to field emitters: the need is for theories of FEVSC in which the BC at the emitter surface is a field value sufficient to sustain the emission (ion or electron emission) of interest.

For electrons, the first theory of this kind was developed by Stern et al [104] in 1929, and a more comprehensive study was reported by Barbour et al [105] in 1953. An overview of then-existing work on 1D models, and a rationalization of 1D theory, was presented by one of us in 2008 [106]. The main results of this 2008 analysis are presented here (but, to avoid clashes, different notation is used).

A convenient language to discuss FEVSC effects is the following. For a given emission voltage, the Laplace field E_Lp is the surface field in the absence of FEVSC and the Poisson field E_Ps is the field when FEVSC is present. A field reduction factor ${\varTheta}_{\textrm{sc}}$ due to space charge is then defined by

$$\begin{equation} {\varTheta}_{\textrm{sc}} = E_{\textrm{Ps}} / E_{\textrm{Lp}} .\end{equation} \tag{ 80 }$$

The role of FEVSC theories is to determine mathematical expressions for ${\varTheta}_{\textrm{sc}}$; with 1D theories the fields involved are interpreted as emitter apex fields. For simplicity in what follows we omit the subscript 'SC'.

A convenient theoretical form is found by defining a constant here denoted by κ_sc and a dimensionless field-dependent parameter here denoted by ξ (but by ζ in [106]), using the equations

$$\begin{equation} \kappa_{\textrm{sc}} \equiv {\epsilon_0}^{-1} (m/2|ne|)^{1/2}, \end{equation} \tag{ 81 }$$

where m is the mass and ne is the charge of the emitted electron or ion, and

$$\begin{equation} \xi \equiv \kappa_{\textrm{sc}} {|V_{\textrm{e}}}|^{1/2} |J_{\textrm{Ps}}| / |{E_{\textrm{Ps}}|}^{\,2} , \end{equation} \tag{ 82 }$$

where $|J_{\textrm{Ps}}|$ is the magnitude of the related LECD, evaluated for field magnitude $|E_{\textrm{Ps}}|$. The LECD can either be an assumed experimental value, or a value predicted by using some specified expression for LECD. The parameter ξ has been called the space-charge strength.

An analysis of FEVSC electrostatics analogous to that made by Barbour et al [105] yields the equation

$$\begin{equation} 9 \xi^{\,2} {{\varTheta}}^{\,2} - 3{{\varTheta}} - 4\xi + 3 = 0, \end{equation} \tag{ 83 }$$

called in [106] the dimensionless planar FEVSC equation.

This has a general solution that can be implemented on a spreadsheet, but the mathematics is messy because there is a branch point at $\xi = 1/2$ (see [106] for details). Several regimes can be identified where simplified equations apply. Those most of interest are as follows.

In the low-ξ regime, the general solution reduces to the Stern et al result (see [104], equation (13))

$$\begin{equation} {\varTheta} = 1 - (4/3) \xi + 3{\xi}^{\,2} - 8{\xi}^3\approx 1 - (4/3) \xi, \end{equation} \tag{ 84 }$$

where the second approximation is valid if ξ is sufficiently small.

In the very-high-ξ regime the general solution reduces to

$$\begin{equation} {\varTheta} = (2/3){\xi}^{-1/2} {[1-3/(4\xi)]}^{1/2}.\end{equation} \tag{ 85 }$$

This very-high-ξ result is the FEVSC strong space-charge limit. The formula has the form [$(2/3){\xi}^{-1/2}$] (see [106]) of the strong space-charge limit for thermally emitted vacuum space-charge (Child's Law), multiplied by the factor in square brackets, which gets closer to unity as ξ increases.

The theory in [106] can be used make an estimate of when FEVSC effects are likely to become significant. From equation (18) in [106] it can be found that ${\varTheta} = 0.90$ corresponds to ξ ≈ 0.08845. Bearing in mind that, in a 1D model

$$\begin{equation} |V_{\textrm{e}}| = |E_{\textrm{Lp}}| d_{\textrm{gap}} = |E_{\textrm{Ps}}| d_{\textrm{gap}} / {\varTheta} , \end{equation} \tag{ 86 }$$

where the separation between the planar electrodes is here denoted by d_gap, equation (82) can be re-arranged into the form

$$\begin{align} |J_{\textrm{Ps}}| / |{F_{\textrm{Ps}}|}^{3/2} & = (\xi {{\varTheta}}^{1/2} \kappa_{\textrm{sc}}^{-1}) /d_{\textrm{gap}}^{1/2} \nonumber\\ & \ {\approx (4.41 \times 10^{-7} \; \textrm{m}^{1/2}}) / {d_{\textrm{gap}}^{1/2}}, \end{align} \tag{ 87 }$$

where we have used the fact that the value ⁴ of κ_sc for an electron is $1.904\times 10^5$ A⁻¹ V$^{3/2}$.

For example for $d_{\textrm{sep}} = 1$ mm, using the 1956 MG FE equation and φ = 4.50 eV, this yields the result that there is a FEVSC-induced 10% reduction in LECD at a field around 5 V nm⁻¹, at which point the LECD is about $5\times 10^{9}$ A m⁻².

Most LAFEs involve irregular arrays of point-like emission sites, but there has been some past interest (e.g. [107]) in arrays of emitters that have the shape of 'walls', 'wedges', 'knife-edges' or 'blades'. We refer to these generally as blade emitters. Practical blade emitters have finite length, and may have rounded ends, but theoretical modelling starts by considering blade emitters of infinite length.

Recently, the growing interest in the possibility of FE from the edges of graphite flakes has been recognised in the publication of a new handbook [14], containing much useful material, in particular an overview of emission from graphene edges by Purcell et al [108].

Emission theory for blade emitters differs from standard forms of emission theory, because the effects of quantum confinement in the direction 'across the blade' need to be taken into account. Early investigations were carried out by Li et al [109, 110]; more recent theoretical developments include [111] and are discussed in [14], in particular in [112].

Early treatments of the electrostatics of single-blade emitters include [113, 114]. Sun and Ang investigated blade emitters with a cusp-like cross-section [115]. Tang et al investigated depolarization effects related to two parallel blade emitters [116]. Filippov et al have investigated straight and rounded blade emitters with several different cross-sectional profiles [117, 118]. There remains scope for further work on the electrostatics of blade emitters, and for a more detailed review in due course.

In principle, line-charge arguments can use either a linear array of point charges (and/or dipoles), or a continuous linear distribution of charge. The aim is to choose a charge distribution such that one of the equipotentials models the emitter-shape of interest. Potentials and fields at other positions are then determined by the chosen charge distribution. If an exact emitter shape is required, choosing the related charge distribution can be non-straightforward.

Line-charge arguments can be applied in both field electron and field ion emission, but have found more use in FE. The theory of applying line-charge distributions to needle-shaped emitters and to field emission has been discussed, for example, by Harris et al [119] and Zheng et al [120].

The first people to have used line-charge arguments in FE seem to have been Vibrans [53, 54] and Egorov and co-workers [121]. They have been used extensively in the work of Harris et al (e.g. [119, 122, 123]), and also by Biswas et al (e.g. [124, 125]).

Holographic-type electron microscopy experiments have been carried out in order to verify the form of line charges, and find [120] that an extra charge near a rod/needle apex is required: this finding is consistent with the basic version of the FSEPP model, where only a single point charge is used to model the emitter.

In general terms, the accuracy of the line-charge models lies between the FSEPP model and the FEM simulations. For a short period, before the recent development of highy accurate FEM simulations described earlier, line-charge models were useful in exploring the two-emitter interaction. Their main usefulness now looks likely to be in niche applications where FEM simulations are difficult to implement. A particular niche may be the theory of small finite arrays of randomly placed emitters. A line-charge approach has also recently been used to model a complicated dielectric structure [126].

At present, most electrostatic modelling in field electron emission is 'above the atomic level', but atomic-level modelling has long been needed in field ion emission to deal with: (a) the interpretation of field ion images of individual atoms (e.g. [127, 128]); (b) the physics of field adsorption (e.g. [129]); (c) the issue of 'which (partially charged) surface atom evaporates next, in the process of FEV (e.g. [34])'; and (d) the locations of the specimen electrical surface (e.g. [130]) and hence of the critical surface inside which ordinary field ionization cannot occur. The electrical-surface location also affects the location of the FE tunnelling barrier [17].

Most of the earlier modelling has been classical in nature, using point charges and dipoles where relevant, but this is now largely being replaced by modern methods based on DFT.

Detailed review of atomistic electrostatic modelling is outside the scope of this article, but we think it useful to indicate the FE contexts in which it is already known that it will be needed in future, as part of the problem of predicting LECDs.

• (a)  
  Understanding the details of FE from CNTs, including the probably important role of patch fields and the possible role of field penetration into the emitter apex. Some particularly relevant discussions are [131–135]. An important conclusion of work so far is the strong indication that better consistency between experiment, classical models and quantum-mechanical models is found if modern DFT calculations are used to derive the induced atomic-scale charge, field and potential distributions, by subtracting off the zero-field results from the finite-field results.
• (b)  
  Understanding the details of FE from flat atomically structured surfaces. There is already some relevant modern work (e.g. [136, 137]), but more is needed. This work needs to be linked with related work on positively charged surfaces (e.g. [138]).
• (c)  
  Investigating the effect of single-atom adsorption and single-atom-defects on otherwise flat surfaces (e.g. [139]).
• (d)  
  Investigating the electrostatics and emission properties of atomically sharp emitters, such as those developed by Fink [140], and of the three-atom tips used in gas field ion sources for scanning ion microscopy (e.g. [141]).
• (e)  
  In both FE and FI emission, investigating sub-atomic level and bond-level emission properties, where these can be now be observed experimentally (e.g. [142–145]).

The exploration of atomic-level and sub-atomic-level FE and FI emission physics seems likely to become a linked growth area in these theoretical specialities. There also seems a need for a 'properly joined up' theory of charged surfaces that relates the classical and quantum mechanical approaches.

A question of basic interest is what is the total field enhancement generated by a two-stage structure, in which one protrusion stands on the apex of another. Experimentally, it is known (e.g. [146–148]) that such arrangements can generate a high total apex FEF. The basic physics can be understood qualitatively by considering PPP geometry and using FSEPP models for two HCP posts.

A 'large' post, of height $h_{\textrm{b}}$ and apex radius $r_{\textrm{b}}$, stands directly on the emitter plate; a 'small' post, of height $h_{\textrm{a}}$ and apex radius $r_{\textrm{a}}$, stands on the apex of the large post. Let $E_{\textrm{B}}$ denote the (inter-)plate field (i.e. the electrostatic field at the emitter plate surface in the absence of both posts). Let $E_{\textrm{b}}$ denote the apex field for the large post, in the absence of the small post, and let the corresponding PFEF $\gamma_{\textrm{Pb}} \equiv E_{\textrm{b}} / E_{\textrm{B}}$.

If the small post were standing directly on the emitter plate (in the absence of the large post), then let its apex field be $E_{\textrm{a}}^{(1)}$ and the related apex FEF be $\gamma_{\textrm{a}}^{(1)} \equiv E_{\textrm{a}}^{(1)} /E_{\textrm{B}}$. When the small post is standing on the apex of the large post, then let its apex field be $E_{\textrm{a}}^{(2)}$ and the related single-stage FEF be $\gamma_{\textrm{a}}^{(2)} \equiv E_{\textrm{a}}^{(2)} /E_{\textrm{b}}$.

By definition, the total apex PFEF for the two-stage structure is $\gamma_{\textrm{Pa,tot}} \equiv E_{\textrm{a}}^{(2)}/E_{\textrm{B}}$. Expansion yields

$$\begin{equation} \gamma_{\textrm{Pa,tot}} \equiv E_{\textrm{a}}^{(2)}/E_{\textrm{B}} = (E_{\textrm{b}}/E_{\textrm{B}}) (E_{\textrm{a}}^{(2)}/E_{\textrm{b}}) = \gamma_{\textrm{Pb}} \times \gamma_{\textrm{a}}^{(2)} .\end{equation} \tag{ 88 }$$

Thus, when the fields E_B, E_b and $E_{\textrm{a}}^{(2)}$, and related apex FEFs, are defined as above, then it follows from these definitions alone that the total apex FEF is found by multiplying the separate FEFs for the two stages. we refer to this as two-stage field enhancement. Analogous multi-stage formulae can also be given.

In the limiting special case where $h_{\textrm{a}}\ll r_{\textrm{b}}$, and the situation can be treated as if the small post were standing on a planar surface, this relation reduces to

$$\begin{equation} \gamma_{\textrm{Pa,tot}} = {(\gamma_{\textrm{Pa,tot}})}^{\textrm{SchC}} \approx \gamma_{\textrm{Pb}} \times \gamma_{\textrm{a}}^{(1)} = \gamma_{\textrm{Pb}} \times \gamma_{\textrm{Pa}} .\end{equation} \tag{ 89 }$$

It is thought that Schottky [149] was the first person to consider the theory of two-stage electrostatic effects of this general kind, and equation (89) is often now called Schottky's Conjecture. Schottky's mathematical analysis was, in fact, somewhat different, and there seems to be a discrepancy in his published equation (see page 78 of [149]).

In reality, Stern et al [104] seem to have been the first in FE to use the physics of formula (89), when they argued that putting a small hemisphere on top of a much larger hemisphere would generate a total apex FEF of 9.

A matter of some interest is how to establish a condition for breakdown of equation (89). The most reliable approach is via finite element methods (e.g. [150, 151]). However, the issues involved can be explored qualitatively by using the idea of differences in 'driving potential difference' (as discussed earlier) to compare the behaviours of a 'small' HCP-model post (of height h_a and apex radius r_a) when: (a) it is located in an uniform applied field; and (b) it stands on the apex of a 'large' HCP model post (of height h_b and apex radius r_b).

In case (a), the driving potential difference between the top and bottom of the small post is $\Delta {{\varPhi}}_1 = -E_{\textrm{b}}h_{\textrm{a}}$, where the post is treated as in an uniform field equal to the field E_b at the apex of the large post.

In case (b) the FSEPP model is used to represent the electrostatic effects of the large post, by placing a charge Q at the sphere centre. The driving potential difference $\Delta {{\varPhi}}_2$ for the small post is then estimated by

$$\begin{align} \Delta {{\varPhi}}_2 & \approx \left[ Q/(4\pi\epsilon_{0}) \right] \left[ 1/(r_{\textrm{b}} + h_{\textrm{a}}) - 1/r_{\textrm{b}} \right] \nonumber\\ & = E_{\textrm{b}} r_{\textrm{b}} \left[ 1/( 1+ h_{\textrm{a}}/r_{\textrm{b}}) - 1\right] \nonumber\\ & \approx E_{\textrm{b}} \left[ - h_{\textrm{a}} + (1/2)(h_{\textrm{a}}^{\,2}/r_{\textrm{b}}) \right].\end{align} \tag{ 90 }$$

Hence the fractional difference N between the cases is

$$\begin{equation} N \equiv (\Delta {{\varPhi}}_2 - \Delta {{\varPhi}}_1) / \Delta {{\varPhi}}_1 \approx -h_{\textrm{a}}/2r_{\textrm{b}}.\end{equation} \tag{ 91 }$$

If, say, we wish the error to be less than 10% ($|N|\lt0.1$), then this model suggests the criterion

$$\begin{equation} h_{\textrm{a}} \lt \; \approx \; 0.2 \;r_{\textrm{b}}.\end{equation} \tag{ 92 }$$

Figure 18 shows the results of a numerical calculation [150] for this two-HCP-stage model. This shows that, in this geometry, significant breakdown of the conjecture in fact occurs for $h_{\textrm{a}} \gt \; \approx \; 0.1 \; r_{\textrm{b}}$. The difference from equation (92) is due to the approximations in FSEPP-type models.

Zoom In Zoom Out Reset image size

Broadly similar conclusions were reached by Harris et al using line-charge models [152].

An equivalent condition for a blade-type emitter was derived by Miller et al [153].

Variants of the analytical theory given above can be developed for cases where the small protrusion stands, not at the apex of the large protrusion, but at some other point on the protrusion.

Various other methods of investigating two-stage field-enhancement effects and/or the limitations of equation (89) have been explored in the literature, for example [113, 146, 148, 151–160].

More generally, it is thought that two-stage field enhancement plays an important role in many physical situations, particularly in some electrical-breakdown situations, and can contribute to the development of highly effective technological field emitters. A particular example of this is the 'looped fibre + fibrils' configuration discussed in [147], where experimental data-analysis that passes an Orthodoxy Test validity check appears to suggest that a two-stage field enhancement factor of order 10⁵ has been achieved.

Collector adjacency effects occur when the counter-electrode (here called the collector) is sufficiently close to the emitter that this affects the values of the apex local field and related field enhancement factors. Adjacency effects have been widely discussed in FE literature (e.g. [55, 78, 148, 161–170]), usually for PPP geometry and HCP model emitters, though also for other models.

In PPP geometry, adjacency effects have been discussed for both the relevant plate-FEF and the relevant gap-FEF. As discussed below, these two FEF-types have qualitatively different behaviour.

It is simplest to consider an electronically ideal emitter of constant height h and constant apex sharpness ratio $\sigma_{\textrm{a}} = h/r_{\textrm{a}}$. One needs to establish the effect of decreasing the gap length d_gap between the emitter apex and the collector, but there are options for discussing this.

Our initial approach is to reduce the interplate separation d_sep whilst keeping the interplate field magnitude $|E_{\textrm{P}}|\,[ = |V_{\textrm{P}}|/d_{\textrm{sep}}]$ constant, by appropriately reducing the interplate voltage magnitude $|V_{\textrm{P}}|$. In these circumstances the magnitude $|E_{\textrm{a}}|$ of the apex local field would be constant, independent of d_sep, if there were no adjacency effects. In practice, as d_sep is reduced from a large value, a point is reached where $|E_{\textrm{a}}|$ begins to increase, and the situation moves towards a narrow-gap regime where $d_{\textrm{gap}} \ll r_{\textrm{a}}$ and

$$\begin{equation} |E_{\textrm{a}}| = |V_{\textrm{G}}|/d_{\textrm{gap}} = |V_{\textrm{P}}|/d_{\textrm{gap}} . \end{equation} \tag{ 93 }$$

Here, $|V_{\textrm{G}}|$ is the magnitude of the voltage between the emitter apex and the adjacent collector; for an electronically ideal system this equals the magnitude $|V_{\textrm{P}}|$ of the interplate voltage.

Since the apex plate-FEF $\gamma_{\textrm{Pa}} = |E_{\textrm{a}}|/|E_{\textrm{P}}|$, it follows that γ_Pa has the same behaviour as that just described for $|E_{\textrm{a}}|$.

If, instead, one holds the interplate voltage constant, then the behaviour of the plate-FEF is the same as just described, but the apex local field magnitude $|E_{\textrm{a}}|$ now has a more complicated dependence on the plate separation d_sep.

In all the above cases, the point at which a 'significant' dependence on d_sep begins to develop (for likely practical values of $h/r_{\textrm{a}}$) is usually stated to be around $d_{\textrm{sep}} \lt h/2$ to $d_{\textrm{sep}} \lt h/3$. However, it is expected that the details will depend somewhat on the shape of the emitter apex, and in particular on the apex radius. If the issue is important, then specific FEM-based calculations are recommended.

In our view, the approach just described, in terms of plate-FEFs, is often preferable. However, there is an alternative approach that starts from the properties of the gap. In the narrow-gap regime, the apex field is given by equation (93). As the gap length d_gap increases, this formula needs (for an electronically ideal emitter) to be replaced by

$$\begin{equation} |E_{\textrm{a}}| = |V_{\textrm{G}}|/\zeta_{\textrm{Ga}} = |V_{\textrm{P}}|/\zeta_{\textrm{Ga}} , \end{equation} \tag{ 94 }$$

where the (apex value ζ_Ga of the) gap-voltage conversion length is defined by this equation.

The magnitude of the (mean) gap-field E_G is defined by

$$\begin{equation} |E_{\textrm{G}}| \equiv |V_{\textrm{G}}|/d_{\textrm{gap}} = |V_{\textrm{P}}|/d_{\textrm{gap}} .\end{equation} \tag{ 95 }$$

Hence the apex value γ_Ga of a gap-field enhancement factor (GFEF) can be defined by

$$\begin{equation} \gamma_{\textrm{Ga}} \equiv |E_{\textrm{a}}|/|E_{\textrm{G}}| = d_{\textrm{gap}}/\zeta_{\textrm{Ga}} .\end{equation} \tag{ 96 }$$

As d_gap increases from a very low value, the gap-FEF transitions from being unity to being effectively equal to the large-gap plate-FEF.

In our view, it may be more straightforward to deal directly with the gap VCL, which transitions from being equal to the gap length to being equal to the usual apex local VCL, but formulation in terms of gap-FEFs is more usual.

Various numerical discussions of the transition process, for both plate-FEFs and gap-FEFs, have been given in the citations earlier, but in our view have not yet reached the stage where a quantitatively definitive summary can be presented.

The general effect reported here, namely transition from a 'narrow-gap regime' to a 'wide-gap regime' as the gap length increases, also occurs in point-plane and point-point electrostatic geometries, but detailed discussion is outside the scope of this review.

When two or more post-like emitters stand on a conducting plane the emitters interact electrostatically. For simplicity, most existing numerical discussions consider only the special case of identical, cylindrically symmetrical, post-like emitters. This review also considers only this special case. The ratio $c/h$ of the emitter separation c to the common height h is then of particular relevance. The term 'emitter separation' refers to the distance between the post axes. The ratio itself is called here the separation ratio.

As already discussed, each post, together with its electrical image in the conducting base-plane, forms a finite electrostatic dipole.

The apex-PFEFs for two or more posts are mutually decreased if the post separation is decreased. Our preferred name for this effect is electrostatic depolarization, since the physics of what is happening is a decrease in the magnitude of the finite electrostatic dipole. Depolarization effects associated with point electrostatic dipoles have long been part of the theory of charged surfaces [129, 171]. In FE literature, the finite-dipole effect is also known as 'screening' and 'shielding'.

In the two-emitter situation, the magnitude of the depolarization can be quantified by the (+ve) fractional reduction ($-\delta$) in the apex-PFEF, as given by

$$\begin{equation} -\delta \equiv \frac{(\gamma_{\textrm{Pa}})_1 - (\gamma_{\textrm{Pa}})_2}{(\gamma_{\textrm{Pa}})_1}, \end{equation} \tag{ 97 }$$

where the subscripts '1' and '2' label the one-emitter and two-emitter situations, respectively. The value of ($-\delta$) depends on the separation c of the emitters.

Depolarization effects also influence field and FEF values at emitter surface positions other than the apex, but the numerics of this has not yet been systematically explored, and such effects are outside the scope of this review.

As discussed by Forbes [52], two separate physical effects contribute to the reduced apex-PFEF value. The primary cause is so-called charge-blunting, which results from the 'driving potential difference' term $(\Delta {\varPhi})_{\textrm{XD}}$ in equation (43). The name 'charge blunting' is used because the effect of $(\Delta {\varPhi})_{\textrm{XD}}$ is to force charge back from the emitter tip onto the conducting plane.

A much smaller contribution is made by the so-called neighbour field effect, in which the field due to the depolarising dipole forms part of the 'external field contribution' E_Xa in equation (44). Both contributions are discussed in [52]. For simplicity, the analytical discussions here mostly disregard the neighbour field effect; however, it is automatically included in the FEM electrostatic simulations.

If ($c/h$) is relatively small (${\approx} 1$, or less), the depolarization related to charge blunting must be described directly in terms of contributions to the driving potential difference. This approach can be used at all separations, as in [52]. But, if the emitters are reasonably well separated, then the effect is conveniently described as due to a depolarizing field E_D. This opposes the inter-plate field and makes, to the driving potential difference, a contribution $-E_{\textrm{D}} h$ that reduces the total magnitude of $(\Delta {\varPhi})_{\textrm{XD}}$.

If the emitters are very well separated, then the theory effectively reduces to the mutual depolarizing effect of point dipoles: this is standard electrostatics and is used (for example) in the classical theories of charged surfaces (e.g. [172]) and of field adsorption, (e.g. [129]), as employed in the context of field ion microscopy.

When more than two emitters are present, the depolarization effects are basically additive. For any particular emitter under analysis, all emitter pairs need to be considered, but with contribution magnitudes that depend on the pair separation. With small clusters or finite arrays, emitters near the centre of the cluster experience greater total depolarization than those near the edges of the cluster or array, due to the geometry of the situation. Thus, field enhancement factors are higher near the edges of the cluster or array.

This effect can be seen in the '7-emitter cluster' figure below, and has also been found in earlier numerical simulations, for example those of Murata et al [173], Bocharov and Eletskii [174] and Smith and Silva [175], and in the line-charge models used by Harris et al (e.g. [123]). Its existence is well established experimentally [private discussion at conferences], though sparsely described in FE literature. This edge effect can be a hindrance in the design of practical devices. Sometimes, steps to mitigate its consequences have been explored (e.g. [123, 176, 177]).

An implication of this theoretical additivity is that the two-emitter interaction is the basic building block for cluster theory and array theory and requires careful attention. This is carried out in sections 10.2–10.4. Depolarization theory relating to larger clusters and arrays of post-like emitters is not yet fully established, and is too large a topic to be included in this review; section 10.6 provides some introductory remarks.

Numerical explorations of the two-emitter situation were initially made by Harris et al [119, 122], using line-charge emitter models, and led to diagrams broadly similar to figure 23. Following earlier empirical thinking (e.g. [178]), they interpreted their work in a manner that implies an exponential fall off in ($- \delta$) with separation, at large separations. They also found [122] that, at very small separations, where the value of c was only a few post radii and was reducing, the apex-PFEF went through a minimum and began to increase again.

The later exploration by Forbes [52], using two FSEPP model emitters, suggested that (apart from the regime at very small separations) there were at least two ($c/h$) regimes. At large ($c/h$) values, ($-\delta$) diminished with the inverse cube of distance. However, as ($c/h$) moved towards unity, there was a change-over to a regime where ($c/h$) changed more slowly. This is illustrated below. A second discussion (arXiv:1803.03167), which was later confirmed in further papers [73, 179], firmly associated the inverse-cubic regime with the electrostatics of finite dipoles.

The two regimes were originally called the 'closely separated' and 'widely separated' regimes. We consider it better to call the first of these regimes (where $c/h$ is roughly between 0.1 and 1) the moderately separated regime, and to suppose the existence of a narrowly separated regime where $c/h \lt {\sim}0.1$ (but $c/r_{\textrm{a}} \gg 1)$.

At very small values of $(c/h)$, de Assis and Dall'Agnol [180] identified a regime where the position of maximum field on the emitter tip moves away from the axis. This regime has been called the close proximity regime and is discussed in section 10.4. This regime presumably either overlaps with or is identical with the narrowly separated regime found with the line-charge models, but details of the relationship have not yet been explored.

In the ideal world one needs, for each of the four separation regimes, physically derived analytical formulae for γ and/or $- \delta$. In practice these formulae will not be accurate, so one also needs both a numerical treatment and a fitting formula. In practice, at present, none of the analytical formulae or fitting formulae cover all four regimes, so one needs to specify ranges of applicability for formulae.

As shown below, a further complication is that values of γ and $- \delta$ will also depend on the apex sharpness ratio, and simulation reports need to take this into account.

All these things will be different in detail for posts of different shapes. We therefore concentrate on the behaviour of the well-studied HCP emitters, on the assumption that behaviours for other post shapes will be qualitatively similar.

However, as things stand at present, even for HCP model emitters our quantitative understanding of the two-post interaction (although well advanced) is not yet complete. Section 10 reports the present situation. As a preliminary, section 10.2 presents a physical explanation of electrostatic depolarization.

Electrostatic depolarization is a physical effect that appears to be unfamiliar to many people. A physical explanation may be helpful.

Consider two identical conducting posts that stand on a grounded horizontal flat conducting plane and are situated in a vertically uniform, macroscopic field E_P. Each post develops a charge-distribution near its tip, and there is a corresponding image distribution on the opposite side of the grounded plane.

In the lowest order of approximation, these two distributions can each be represented by a finite electrostatic dipole of moment p. Thus, from the electrostatic point of view, we are considering two parallel dipoles that tend to mutually depolarize each other, as illustrated in figure 19. In the lowest order of approximation (which is adequately valid if $(c/h)$ is sufficiently large, each dipole generates, at the base of the other post, a depolarizing field E_D given by:

$$\begin{equation} E_{\textrm{D}} = - \frac{p}{4\pi \epsilon_0 c^3}.\end{equation} \tag{ 98 }$$

There is an issue of what are the best units to use for measuring p. This is addressed in the appendix.

Zoom In Zoom Out Reset image size

If c is sufficiently large, then the depolarizing field strength varies relatively little along the length of the post, and we can define an effective electrostatic field E_eff acting on the post by

$$\begin{equation} E_{\textrm{eff}} = E_{\textrm{P}} + E_{\textrm{D}}.\end{equation} \tag{ 99 }$$

The polarization of each emitter is actually induced by E_eff rather than E_P. Hence,

$$\begin{equation} p = \alpha E_{\textrm{eff}}, \end{equation} \tag{ 100 }$$

where α is the effective polarizability of each emitter. As a result, combining equations (98)–(100) yields

$$\begin{equation} E_{\textrm{eff}} = {{\varTheta}}_{\textrm{D}} E_{\textrm{P}} \equiv \frac{E_{\textrm{P}}}{1 + \alpha/4\pi \epsilon_0 c^3}, \end{equation} \tag{ 101 }$$

where the correction factor ${\varTheta_{\textrm{D}}}$ is defined via equation (101).

In the theory in section 5.3, E_B there needs to be identified as the plate-field E_P here. Thus, in the context of the $(\Delta {\varPhi})_{\textrm{XD}}$ term in section 5.3, $-E_{\textrm{B}} h$ there is replaced by the expression $-E_{\textrm{eff}} h$ [$= - {{\varTheta}}_{\textrm{D}} E_{\textrm{P}} h]$, and the outcome is the following formula for the fractional reduction $(- \delta)$ in apex-FEF:

$$\begin{equation} (-\delta) = 1 - {{\varTheta}}_{\textrm{D}} = \frac{\alpha}{\alpha + 4\pi \epsilon_0 c^3}.\end{equation} \tag{ 102 }$$

In the limit of large separation, this becomes an inverse-cubic dependence on separation. This prediction has been verified by FEM simulations, as illustrated in following discussion.

Clearly, result (102) is not strictly valid if the depolarizing field varies along the post. In such circumstances, it is better to treat the depolarizing dipole as a pair of negative and positive charges, as is done in the Forbes FSEPP-model approach [52], and directly develop an expression for the driving potential difference $(\Delta {\varPhi})_{\textrm{XD}}$ in equation (43).

This FSEPP-model approach yields a general expression (for the charge-blunting contribution) of the form

$$\begin{equation} (- \delta) = \frac{C_1}{C_1 + C_0} \approx C_1 , \end{equation} \tag{ 103 }$$

where C₀ and C₁ are geometrical parameters defined algebraically in [52]. In the limit of large separation (i.e. in the 'widely separated regime') $-\delta$ is given adequately by

$$\begin{equation} (- \delta) \approx 2r_{\textrm{a}} {\ell}^{\,2} / c^3\approx 2\times \left( \frac{h}{r_{\textrm{a}}} \right)^{-1} \times \left(\frac{c}{h} \right)^{-3} , \end{equation} \tag{ 104 }$$

where r_a denotes the radius of the floating sphere, and $\ell = h-r_{\textrm{a}}$ is the height of the centre of the floating sphere above the ground plane. It can be shown that this expression is exactly equivalent to $\alpha/4\pi \epsilon_0 c^3$.

At somewhat smaller separations, but in conditions where $c \gg r_{\textrm{a}}$ (i.e. in the 'moderately separated regime') $-\delta$ is given adequately [52] by

$$\begin{equation} (- \delta) \approx r_{\textrm{a}} / c \approx \left( \frac{h}{r_{\textrm{a}}} \right)^{-1} \times \left(\frac{c}{h} \right)^{-1}. \end{equation} \tag{ 105 }$$

These regimes, together with a plot of the exact variation of C₁ are shown as a ln-ln plot in figure 20, for an illustrative emitter with $\ell/r_{\textrm{a}} = 100$. The analytical models display the underlying physics of the two-emitter situation, but for accurate numerics it is best to use FEM electrostatic simulations. Figure 21 shows FEM-simulation results for geometries that include that considered in figure 20. The inverse-cubic fall-off in ($- \delta$) at large $c/h$ values is well confirmed, but the 'moderately separated' regime is poorly defined and has a limited range, and in the 'narrowly separated regime' $( - \delta)$ goes through a maximum as $c/h$ decreases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Equations (104) and (105) predict the physical effect that $(-\delta)$ will depend on the apex sharpness ratio $h/r_{\textrm{a}}$, in both the 'moderately separated' and 'widely separated' regimes. Since this applies to the two-emitter situation, it will also apply to clusters and arrays. However (as noted in our earlier paper [73]), this factor (or dependence) does not clearly appear in most earlier work on electrostatic depolarization, either in analytical discussions of multi-emitter situations or in line-charge models of such situations.

In the circumstances of figure 21 the effect is relatively small. For HCP-model emitter and its ring-like image with respect to the rotated right-hand boundary figure 22 shows that it appears more obviously in the widely-separately regime when a greater range of values of $h/r_{\textrm{a}}$ is considered.

Zoom In Zoom Out Reset image size

The implication of this confirmation is that all earlier analytical treatments of multi-emitter interactions (except, perhaps, those in [52]), should be treated with moderate reservations, and that all earlier line-charge treatments and FEM simulations should be treated as applying only to the apex sharpness ratio assumed in the research in question.

In figure 21, the maxima in simulated values of ln$(- \delta)$ imply that the reduced apex-PFEF value goes through a minimum and then begin to increase again, as $c/h$ decreases. As already noted, the existence of this minimum was originally shown by Harris et al [122], using line-charge modelling. They also found, as would be expected, that for an 'internal' emitter in a cluster, the depth of the minimum depends on the number of emitters in the cluster.

Careful analysis [180] has shown that (in the two-emitter situation) a second effect also occurs. This is that the position 'max' of maximum field enhancement moves away from the emitter apex, in a direction away from the other emitter. This effect is demonstrated in figure 23, for a illustrative apex sharpness ratio of 300, and (as already noted) has been called the close proximity effect.

Zoom In Zoom Out Reset image size

This figure also shows that this effect occurs with closely spaced clusters of emitters, with the displacement of the maximum PFEF γ_Pmax always being outwards from the centre of the cluster, and with (obviously) the extent of the reduction in γ_Pmax increasing with the number of emitters in the cluster. With large clusters or large finite arrays of closely spaced emitters, this 'displacement outwards' effect is expected to mainly affect the emitters at the edges of the cluster/array.

The magnitudes of all these effects are significantly influenced by the apex sharpness ratio, although qualitative features remain the same. For details see [180].

Ideally we need to have fitting formulae for all the FEM-simulations results, but this issue remains to be systematically addressed. For the widely separated regime, a formula has been suggested [73] that we write here in the form

$$\begin{equation} (-\delta)_{\textrm{fit}} = \frac{K_1}{K_2 + (c/h)^3}.\end{equation} \tag{ 106 }$$

For two HCP-model emitters, with $h/r_{\textrm{a}} = 51$, figure 24 illustrates the use of this formula. In the range considered ($c/h \gt 1.5)$ the fit appears to be good.

Zoom In Zoom Out Reset image size

The values of K₁ and K₂ will depend on the apex sharpness ratio, but the good fitting properties of equation (106) in the 'widely separated' depolarization regime have now been verified for many different physical systems.

For the hemisphere-on-plane (HSP) emitter model, the best fit is given by $K_1 = (0.9433\pm 0.0006)$, $K_2 = (12.00 \pm 0.02)$. This result has relevance to discussion of the effect of 'bumps' on a macroscopically flat emitter surface.

In principle, four types of multi-emitter situations are of theoretical interest: larger regularly spaced clusters (both linear and non-linear), finite regular arrays, infinite regular arrays, and random clusters and arrays. Although there is a quite a large body of work on regularly spaced finite clusters and arrays, much of it is in need of reassessment, and we have felt it imprudent to provide extensive citations here. An overview of some of this work has been provided by Harris et al [183] and a useful FSEPP-model treatment has been given by Zhbanov et al [184]. What is really needed is a separate critical review that builds on the detailed discussion of the two-emitter situation in this review, and on older work involving arrays of point dipoles, and uses derived insights to re-assess and extend the existing FE 'cluster' and 'array' literature.

Exceptions to this thinking are a few papers (some already cited) that have carried out numerical simulations on finite regular arrays, and some papers on random arrays that explore some of the specific issues involved. For the random arrays, we note in particular the numerical analysis of RB [75], the analysis of Biswas and Rudra [124] and the conceptual discussion by Bieker et al [185].

In this review, which we think is the first of its kind, we have employed systematic terminology and definitions and have described the basic principles and formulae of field emitter electrostatics. The theory has been formulated in a 'polarity independent' manner that applies to both field electron and field ion emitters.

The distinction between electronically ideal and electronically non-ideal emitters has been noted, as has the fact that the 'electrostatics' of electronically ideal emitters needs to be seen as part of the theories of the chemical thermodynamics and statistical mechanics of electrons.

For ideal emitters, when the emitter is in static electrical equilibrium, the Fermi level is constant throughout the emitter. Making the simplifying assumption that the local work function is uniform across all surfaces in the system leads to the requirement that the classical electrostatic potential is constant across (and 'just outside') the surface of the emitter and its support electrode. This constant-potential surface provides the BC both for classical analytical analyses and for electrostatic simulations using finite elements methods.

In this review we have concentrated on two application areas: (a) the electrostatics of the needle-type emitters used in the projection technologies (field electron and ion sources, field electron and ion microscopies and spectroscopies, and atom-probe microscopy); and (b) the electrostatics of the post-like emitters used in basic models of the arrays of emitters used in LAFEs.

In view of the continuing research interest in FE from CNTs and related structures, we have paid special attention to the electrostatics of the HCP emitter model and to the interaction between two posts of this type, and have identified four separation regimes. We have provided a careful review of our 'minimum domain size' finite-element procedure for finding highly accurate values for plate-type field enhancement factors; this procedure is applicable both to the HCP model and to other post-type emitter models.

We have also provided brief accounts of, and some key citations for, some other effects that can influence FEF values.

We also note that we expect the 'conducting post formula' to find many applications in electrostatics in general.

The writing of this review has drawn our attention to aspects of field emitter electrostatics where further basic research and/or a mini-review of existing knowledge could be useful. We think the following topics would be of current interest.

• (a)  
  A more careful treatment and discussion of the correction factor c_a in the conducting post formula, equation (54).
• (b)  
  Review and extension of knowledge about how FEFs vary with location on the emitter surface, particularly near the emitter apex.
• (c)  
  Review of two-stage field enhancement and further FEM investigations into the limits of validity of Schottky's conjecture.
• (d)  
  Further detailed analysis of the two-identical-post situation: 1. to explore the narrowly-spaced and close-proximity separation regimes more carefully; 2. to establish the influence of the apex sharpness ratio; and 3., if possible, to establish a fitting formula or formulae that cover all regimes.
• (e)  
  Systematic treatment and review of the theory of arrays, initially infinite and finite regular arrays, later random arrays and clusters.
• (f)  
  More detailed analysis and review of collector adjacency effects.
• (g)  
  Exploration of depolarization effects amongst emitters with different apex sharpness ratios.
• (h)  
  Review of the theory of LAFEs, particularly when statistical effects are taken into account that result from a distribution of emitter geometrical parameters.
• (i)  
  A review of the electrostatics of blade geometries.
• (j)  
  Review of the electrostatics of point-plane situations (including 'ball-probe' geometries).
• (k)  
  Development of the electrostatics of three-electrode devices, such as nanoscale vacuum channel transistors.
• (l)  
  More careful exploration of the effects of variations in local work function, and of resulting patch fields.

All the above need to be done for metal systems. There is also a need for equivalent results where non-metals are involved, in particular for the electrostatics of semiconductor posts, and for point-plane analyses where the plane is the front surface of a semiconductor.

All the above involve the theory of classical conductors and semiconductors. As already noted, there is also a need for atomistic-level electrostatics. Much attention has already been given to CNTs and some to other carbon-related materials, but further research is needed, in order to make improved connection between FE theory and experiment. The electrostatics of atomically structured metal surfaces also needs more detailed research, but probably the immediate need is for a brief review that brings together existing results from both field electron and field ion emission.

A particular problem with DFT calculations on CNTs is the limited number of atoms (and hence the limited emitter size) that can be used. In this context we draw attention to the promising hybrid methodology introduced by Edgcombe et al [186–188], in which classical FEM modelling is used to create initial BCs for an atomistic model representing the near-apex part of a CNT of physically realistic dimensions.

Beyond this there is likely to develop a need for sub-atomic-level electrostatics, in order to help interpret related experiments, but this is largely unexplored research territory.

In short, there is a vast amount of future research that needs to be done in the context of field electron and ion emitter electrostatics. This review has been able to cover only a limited part of the topic. But our hope is that our work can provide both a stimulation and a partial basis for future field emitter electrostatics research. We also hope that our work can contribute more generally to electrostatics knowledge.