General scaling laws of space charge effects in field emission

The characteristics of field electron and ion emission change when the space charge formed by the emitted charge is sufficient to suppress the extracting electric field. This phenomenon is well described for planar emitting diodes by the one dimensional (1D) theory. Here we generalize for any 3D geometry by deriving the scaling laws describing the field suppression in the weak space charge regime. We propose a novel corrected equivalent planar diode model, which describes the space charge effects for any geometry in terms of the 1D theory, utilizing a correction factor that adjusts the diode's scaling characteristics. We then develop a computational method, based on the Particle-In-Cell technique, which solves numerically the space charge problem. We validate our theory by comparing it to both our numerical calculations and existing experimental data, either of which can be used to obtain the geometrical correction factor of the corrected equivalent planar diode model.


I. INTRODUCTION
The current extracted from an electron emitting cathode or an ion emitting anode can be increased by applying higher electric fields, increasing the emitter temperature or irradiating it with light. However, the current density cannot increase beyond a certain limit, due to the space charge effect [1,2], i.e. the suppression of the extracting field due to the space charge formed by the emitted particles. The space charge (SC) effect plays a very significant role in all forms of electron and ion sources [3][4][5][6][7][8][9] and is of paramount importance for the understanding of the ignition of vacuum arcs [10][11][12][13].
Despite this importance, there is no general three dimensional (3D) theory describing the SC effect for curved emitter geometries, mainly due to the high complexity of the problem. In order to calculate the SC effects one has to solve self-consistently three coupled problems: the Poisson equation, the continuity equation, and the electron emission equation. This is possible by utilizing various numerical methods, such as Particle-In-Cell (PIC) [14][15][16], molecular dynamics [17], or the point charge method [18].
Analytical solutions of the SC problem exist only for 1D geometries, where the continuity equation has a simple solution. For instance, Child [1] and Langmuir [2] solved the problem for a planar diode geometry and the special case of fully SC-limited charge flow, i.e. with a boundary condition of zero electric field at the emitting electrode. To distinguish this case from the general problem with a non-zero field at the emitter, we shall call it special SC problem. Later on, Langmuir and Blodgett developed analytical solutions of the special SC problem in the form of a series for non-planar -but, still 1D-geometries of concentric cylindrical [19] and spherical [20] * akyritsos1@gmail.com electrodes. In a recent work [21], it was shown that both problems yield general scaling laws similar to the planar diode.
The general SC problem, which is relevant for field electron and ion emission, was solved for the planar case by Stern et. al. [22]. This solution has been widely used as a reference model to estimate the SC effects on field emission [5,[23][24][25], due to its simplicity and intuitiveness. The connection of the planar diode model to an experimental emitter geometry is typically done via the equivalency introduced by Barbour et. al. [8], i.e. by assigning the real values of voltage and cathode field to a planar diode.
However, PIC simulations [13,15,16] have clearly shown that this equivalence is not always valid, as it tends to significantly overestimate the SC suppression of the field, and thus may lead to wrong values of emission current from a 3D emitter. In this work, we tackle this problem by developing a general 3D theory for space charge suppressed emission. In section II A we derive the general scaling laws for the emission behavior in the weak SC regime, by introducing the corrected equivalent planar diode (CEPD) model and showing that any 3D emitter is equivalent, regarding SC, to a planar diode of certain characteristics, determined by a single geometrydependent correction factor ω. In section II B we derive ω for the spherical and cylindrical diode problems. In section III A we describe a general numerical method to obtain ω for any 3D electrode geometry. We finally validate our theory by comparing to both numerical calculations and existing experimental data in section IV, showing that the correction factor can be obtained for a certain geometry by fitting to either numerical calculations or experimental data.

A. General formulation
The standard formulation of the space charge problem adopted already since Child's work [1], assumes a continuous charge density distribution of the emitted charge ρ(r) and zero initial velocity for the emitted charged particles. Some studies have considered the case of non-zero initial kinetic energy [18], but since the latter is of the order of a few eV, we shall consider it negligible for the cases of emission under high electric field. Under these assumptions, the Poisson equation becomes with boundary conditions Φ = 0 at the emitter and Φ = V at the collector electrodes. Here Φ is the electrostatic potential, k = −1 0 m/2q is a constant that depends on the mass to charge ratio m/q of the emitted particles and the dielectric constant 0 ; V is the applied voltage. The current density distribution J(r) obeys the continuity equation ∇ · J = 0, with boundary condition J(r s ) = J s (r s )n(r s ), where J s (r s ) is the emitted current density andn(r s ) is the normal unit vector at the emitting surface point r s . In order to obtain J and Φ, the above equations have to be solved self-consistently, along with the surface emission laws that give J s (r s ) as a function of the local electric field F (r s ).
The above problem cannot be solved analytically, apart from the planar case, where J is constant. In this case, the solution of the 1D Poisson equation for a planar emitting diode with a gap distance d, voltage V , current density J and cathode field F , yields [22] 6(kJ) 2 In case the current density J supplied at the cathode depends on the cathode field F (e.g. field electron or ion emission), eq. (2) needs to be solved self-consistently with the law that gives J(F ), e.g. the Fowler-Nordheim (FN) law [26,27] or its modern generalizations [28][29][30][31]. Barbour et. al. [8] used approximate methods to achieve this self-consistency, while modern iterative numerical methods allow to find a more accurate solution at relatively low computational costs [12]. Eq. (2) can be simplified by introducing the reduced dimensionless variables [23][24][25] In (3), θ is the "field reduction factor", i.e. the factor by which the field has been reduced from the "Laplace field" F L ≡ V /d due to the SC. ζ is indicative of the "space charge strength", since it is evident from (3) that θ reduces from 1 to 0 (F from F L to 0) as ζ increases from 0 to 4/9, where the Child law [1] limit occurs. Eq. (3) can be solved analytically [24]; yet, it is more convenient to express the physical solution in a perturbation series around ζ = 0, yielding Truncating this series at the linear term yields a good approximation with an error of less than 10% for ζ < 0.25 and θ > 0.6. This "weak SC regime" covers most practical cases in field electron and ion emission, where typically J s 10 12 A/m 2 ) [4,6]. The generalization of this scaling behavior of θ(ζ) for non-planar geometries is the main purpose of this article.
Consider an emitter with a general 3D geometry and a point of interest r s at the emission surface, with local field F . We shall express the Poisson equation in terms of the reduced variables φ ≡ Φ/V , J(r) = J s ξ(r), where J s = J(r s ) and ξ(r) is a unitless variable expressing the variation of the current density from r s . Assuming that the distribution of the surface emission J s does not vary significantly and using the linearity of the continuity equation, we can approximate that ξ(r) depends only on the emitter geometry and not on J s . Finally, we use the reduced space coordinater = r/χ, where χ is the "conversion length" χ ≡ V /F L (F L is the Laplace (J = 0) field at r s ). Equation (1) becomes where∇ denotes derivatives with respect tor. Note that the reduced positionr is scale invariant, i.e. it remains unchanged under a geometrical scaling r → ar, for any scaling factor a, as F L scales as 1/a and χ as a.
The significance of the above equation becomes evident in view of the correspondence of the parameter in the parenthesis ζ ≡ kJ s χ 2 V −3/2 = kJ s V 1/2 /F 2 L with the space charge strength of the planar diode. The classical planar diode equivalence [8,25] corresponds to inserting the above ζ to eq. (4), or equivalently substituting d with χ in eq. (2). Now we shall derive an asymptotic expansion similar to eq. (4) for a general 3D geometry, showing that the above equivalence is not valid, and propose a new one that holds up to the first order on ζ 1. We express eq. (5) in its integral form, utilizing the Green's function for Dirichlet problems [32] G(r,r ), i.e. the solution of the boundary value problem (BVP) ∇ 2 r G(r,r ) = δ(r −r ) on Ω, G(r,r ) = 0 on ∂Ω, with δ(·) being Dirac's functional and Ω the vacuum domain. Then the solution can be written in the form of a Fredholm integral equation as where φ 0 (r) is the Laplace solution (ζ = 0). For ζ 1, φ(r) can be expanded in an asymptotic series using the Adomian decomposition method [33], as By taking the gradient we obtain a similar expansion for the field reduction factor θ = F/F L The above equation demonstrates the inadequacy of the standard planar diode equivalence of Refs [8,25]. The scaling of θ(ζ) for ζ 1 is the same as the planar diode only if ω = 1, which does not hold in general. However, ω can be incorporated as a correction to the planar model. Thus, for a given geometry and a surface point r s , we define the corrected equivalent planar diode model, for which the SC strength is Practically, the CEPD model can be used by replacing ζ = ζ c in eq. (3), or using the CEPD gap distance d c = ω 2 χ and voltage V c = ω 2 V in eq. (2) and solving for F . Alternatively, one can substitute directly V and d = χ in (2), as in Refs [8,25], but use a corrected current density value J c = ωJ s . An important property of the geometry correction factor ω is its scale invariance, resulting from the fact that all variables in eq. (8) are invariant to geometrical scaling. In other words, ω does not depend on the absolute size of the electrodes, but only on their relative shape, thus simplifying its calculation and tabulation for various geometries.
The usage of a correction factor such as ω has been proposed by Forbes [25]. Here we define and calculate it rigorously. ω can be formally calculated for a given emitter geometry and point r s using (8), which is though quite cumbersome. It is practically more convenient to calculate θ(ζ) numerically (e.g. by PIC) for a given geometry, and fit ω to the results, thus yielding a more general and computationally cheap approximate solution of the SC problem. Before advancing to such simulations, we shall focus on two diode geometries, for which the CEPD can be obtained analytically, giving useful physical insights.

B. Analytical solutions for spherical and cylindrical geometries
The geometries of the spherical and cylindrical diodes, for which the electrodes are two concentric spheres or cylinders correspondingly, have attracted theoretical interest since the early SC studies [19,20,34,35]. Langmuir and Blodgett [19,20] solved the specific SC problem for both geometries. Furthermore, for the spherical case, which has been recently used to model sharp field emitters [5], Aizenberg [35] obtained an approximate solution for the general SC problem. Here we apply the CEPD model and calculate ω for both the spherical and cylindrical geometries, thus obtaining the corresponding scaling laws for the weak SC regime.
The resulting expressions for the CEPD correction factor ω are , for the spherical and cylindrical diodes correspondingly, wherer is the ratio between the emitter and the collector radii and D(·) denotes the Dawson function [36]. The derivation is given in appendix A. It is worth noting that forr → 1, i.e. when the emitter and the collector have similar sizes, ω (S) ≈ 1 + 2 5 (r − 1), and ω (C) ≈ 1 + 1 5 (r − 1). This is expected due to the fact that the spherical and cylindrical diodes become similar to the planar one, when the electrodes have similar radii and the gap distance diminishes. Moreover, forr → ∞, i.e. when the collector electrode becomes infinitely larger than the emitter, ω scales as log(4r) for the spherical case and as 3r(logr) −2 for the cylindrical one. This means that the cylindrical diode deviates from the planar one much faster for increasingr. In order to validate the above results, we solved the Poisson equation numerically by the Runge-Kutta method [37]. The results are shown in figure 1 where we plot the field reduction factor θ vs both the SC strength ζ (dashed lines) and the corrected one ζ c = ωζ (solid lines), for various values ofr, both for the cylindrical and the spherical diodes. It is evident from the dashed θ(ζ) curves, that both the cylindrical and spherical geometries deviate significantly from the planar one and behave very differently depending on the radii ratior. However, after the correction ζ c = ωζ, they all collapse in a curve that is very close to the one of the planar diode (blue). Hence, the CEPD model can describe accurately these two curved geometries for a very wide range of parameters. We note that the field suppression given by the CEPD model starts deviating significantly only for the strong SC regime, θ 0.2, while in the derivation of the model we aimed to fit only the linear term of eq. (8) (dashed line). Such an excellent agreement demonstrates the validity of the model far beyond the weak SC regime, which is limited to θ 0.6. Finally, note that the deviation increases withr, with the increase being much faster for the cylindrical diode.

A. General computational method
In this section we shall generalize the calculation of ω, using a numerical method applicable to any geometry. Our technique is based on the finite element code FEMOCS [38], which has been recently enhanced with PIC capabilities [13]. In this work we do not aim to reproduce the temporal evolution of the SC distribution, but only the steady state, which is reached in a subpicosecond timescale [13]. In the steady state, the electric field distribution is constant and all particles emitted from a given point follow the same path to the collector. Therefore, there is no need to solve the Poisson equation concurrently with the particle movement, but rather obtain a self-consistent solution of the electric field distribution and the particle trajectories. A similar method has been used before by Zhu and Ang [16].
Our method starts by considering no SC and solving the Laplace equation in the vacuum region, using the Finite Element Method (FEM), with a Dirichlet boundary condition Φ = 0 at the emitter and Φ = V at the collector (see figure 2). Then the electric field distribution is calculated on the emitter surface and the field emitted current density is obtained utilizing the electron emission computational tool GETELEC [31], which evaluates the appropriate emission formulas. Then charged superparticles (SPs) are injected from the surface into the vacuum. Unlike our previous work [13], here we inject one SP at the center of each surface element and adjust the corresponding SP weight according to the locally emitted current. Thus the weight of the emitted SP is w sp = J f A f ∆t/e, where J f is the emission current density on a certain face element, A f is its area, ∆t is the path integration timestep, and e the elementary charge.
After injection, the electron paths are followed by numerically integrating Newton's equations of motion, until all electron SPs reach the anode boundary, where they are removed. As in Ref. [13], we use the explicit leapfrog method integration scheme [39]. However, since here we are integrating only the electron paths without a con-current field and emission calculation, we can utilize an adaptive integration timestep. The latter starts at 0.1 fs upon particle injection, and is increased or decreased by 15% (value found to ensure numerical stability) at some of the timesteps, in order to maintain an average of 2 timesteps that the particles stay in the same FEM cell. This adaptive timestep technique ensures that near the emission surface, where the paths have high curvature and the charge density is high, the timestep is sufficiently low to give good accuracy. On the other hand, far from the emitter where the densities are very low and the paths have large radii of curvature, the integration remains at feasible CPU times, even for diode geometries where the electrodes have vastly different length scales (note that in figure 2, R/r 0 2 × 10 5 ). The decrease in the necessary CPU time is several orders of magnitude for such diodes.
At each integration timestep, the SPs give a contribution to the charge density (see eq. (5) and (6) of Ref. [13] for details), thus building up the SC distribution which is inserted in the assembly of the right hand side of the finite element Poisson equation. The latter is then solved obtaining a new electric field distribution. This process is repeated as a fixed point iteration, until convergence is reached, after typically 10-20 iterations. After convergence, the electrostatic potential, field, charge density and current density distributions have been obtained. As a convergence criterion, we demanded that the relative root mean square change on the charge density distribution is less than 10 −4 . We note that for the calculation of an I − V curve, the voltage is gradually increased with small steps, utilizing the converged charge density of the previous step as the initial guess of the next. This ensures the stable and fast convergence of the fixed point iteration method, as the initial guess is always relatively close to the desired convergence.
Finally, in order to obtain the CEPD parameter ω for a given point on the emitter surface, we calculate the distribution of J s , F, and F L on the emitting surface for various applied voltages V , using the numerical method described above. For each point we can then calculate the corresponding J s −F L curves of the CEPD model, by solving eq. (3) self-consistently with the emission laws, using d = χω 2 . Then we optimize the value of ω in order to minimize the error between the CEPD and the numerical curves.

B. Geometrical model for comparison to experiment
Although our theory and computational methods are general with respect to the kind of emitted particles (electrons or ions) and the electrode geometry, we shall now focus on a specific example problem of field electron emission, for which experimental data are available for comparison. In the experiments of Ref. [8], an electrochemically etched tungsten cathode was utilized to take I − V measurements at both low and high currents. The cath-ode was then coated with Ba in various coverages, in order to reduce the work function W and reach SC-relevant current densities at achievable applied voltages.
The shape of the cathodes used by the Linfield group [4,8,10,40,41], including the one of Ref. [8], are generally described by the general sphere-on-a-cone (SOC) model, developed by the same group [42]. In this model, the electrodes are shaped as equipotential surfaces of the electrostatic potential produced by an isolated charged sphere-on-a-cone electrode. Such surfaces are defined by r n − α 2n+1 r −n−1 P n (cos(θ p )) = C, where (r, θ p ) are the spherical coordinates (radius and polar angle correspondingly), P n (·) denotes the Legendre function of the first kind [43] of order n ∈ (0, 1), α is the radius of the sphere on the cone and C is a constant parameter that determines which equipotential surface is defined by the equation. The value of n is determined by the aperture angle of the cone γ, via the relationship P n (cos(π − γ)) = 0.
The shapes of the two electrodes are determined by choosing two values of the parameter C. The latter is determined from the desired radii of curvature at the apex (θ p = 0) of each electrode, r = r 0 for the emitter and r = R for the collector, by evaluating eq. (10). Figure 2 gives a comprehensive schematic of the model, along with the quadrangular tesselation we used to solve the Poisson equation by FEM. The blue line corresponds to the emitter (cathode) surface, while the green one to the collector (anode). In the inset at the upper right corner we zoom (several orders of magnitude) in the emitter apex region, where the virtual SOC electrode that defines the model geometry is shown in a magenta line.
Barbour et. al. [8] took micrographs of their tip and fitted a SOC model to it, however they did not explicitly report the parameters they extracted. For this reason, to simulate their tip, we used the standard parameters described in Ref. [42] n = 0.1 (corresponds to γ = 0.78 o ), α/r 0 = 0.235, R = 6.5 cm, and h = R, apart from the scale parameter r 0 , which was fitted to the experimental I − V data, resulting in a value r 0 = 315 nm. Although a direct comparison of these parameters to the emitter shape extracted by Barbour et. al. from their micrograph is not possible, an indirect comparison gives very good agreement. Our simulated geometry produces a theoretical field conversion factor (eq. (3) in Ref. [42]) β ≡ χ −1 = 0.3617 µm −1 , which is very close to the value 0.37 ± 0.06 µm −1 reported by Barbour et. al. for the geometry they extracted from micrographs. Note that the theoretical value of β is slightly lower than the one calculated numerically, because the theoretical electrode SOC shapes extend to infinity. Nevertheless, to approximate the real electrode shapes, we used a finite height h = R, which produces a slightly higher field enhancement. The value of h was chosen in view of the realistic electrode setup, as depicted in figure 2 of Ref. [8].

IV. RESULTS
Given the geometrical model described in section III B, we used the computational method described in section III A to solve the space charge problem and calculate the emission current density, charge density, and potential distributions, for various values of the applied voltage V . In figure 3 we plot the emitted current density J s and the surface electric field F as a function of the Laplace field F L at a point on the emitter surface with polar angle θ p = 45 o . The numerical results are shown in dots, while the CEPD and EPD calculations are shown in solid and dashed lines correspondingly. We see that the CEPD model follows the PIC results very accurately, while the EPD model deviates significantly. The value of the CEPD correction factor ω, which is calculated by fitting to the PIC results as described in section III A, varies on the emitting surface, mainly due to the variation of the current density distribution ξ(r). However, this variation of ω is small in the emission area (less than 10%), allowing to use a single effective value ω eff that describes the whole emitter geometry, without significant loss of accuracy. The value of ω eff is fitted by minimizing the deviation between the total current -voltage (I − V ) curve as calculated numerically and as estimated by the CEPD model. The total current I is calculated by numerically integrating J s over the emission surface for both cases. In figure 4, we compare the experimental data of Ref. [8] (dots) with our calculations of the total current using the CEPD model with a single effective correction factor ω eff . We chose r 0 = 315 nm by fitting the theoretical curve for work function W = 4.5 eV to the experimental data for a clean tungsten cathode, in the low-field regime where SC effects are negligible. r 0 determines both the conversion length χ (equivalently the enhancement factor β), i.e. the slope of the curve, and the effective emission area, i.e. the vertical shift of the curve. We see that the value r 0 = 315 nm, which yields a calculated conversion length χ = 2.353 µm, produces an almost perfect match to the measurements. Furthermore it is compatible with the shape extracted by the emitter micrographs in Ref. [8]. Although the work function varies even on the surface of a clean emitter, we assumed a uniform effective work function W for all the curves, indicated in the legend. For the clean surface we assumed the standard tungsten value of 4.5 eV, in line with Ref. [8]. The corresponding values for the coated cathodes, for which no prior knowledge is available, were fitted to match the measurements in the low field regime, after choosing r 0 from the clean-surface curve. Finally, we note that for the lowest work function case, the current was multiplied by a fitted correction factor of 0.78, to account for the reduction of the effec-tive emission area due to the increased non-uniformity of the work function, evident from the corresponding micrograph of figure 5 in Ref. [8]. We then used the numerical method described above to calculate the CEPD model correction factor for this emitter geometry, yielding a value ω eff = 0.8. Note that this value is slightly higher than the specific value of ω = 0.75 extracted for the point with θ p = 45 o in figure 3, because ω eff is an effective average of all the surface points.
We see that the CEPD model predicts very accurately the curvature of the experimental plots caused by the SC effects at high fields. Surprisingly, the standard EPD model as introduced in [8], gives also good agreement with the measurements. This is due to the fact that ω eff is very close to unity, resulting in a minor deviation between the EPD and CEPD models, as seen by the hardly distinguishable dashed curves.
FIG. 5. Comparison of I − V curves for emitters of various cone apertures γ as calculated by PIC (dots) and by the CEPD (solid lines) and EPD (dashed lines) models. The calculated correction factor ω eff of the CEPD model is given in the legend for each geometry. The experimental data [8] for the clean emitter are given for comparison.
However, this is rather a coincidence, specific to this particular geometry. In figure 5 we plot I − V curves for different geometries, varying the cone aperture γ, while all other geometrical parameters are kept equal to the ones of figure 4. We see that as γ increases, ω eff decreases along with the field enhancement factor. This results in a significant deviation of the EPD model from the PIC calculations, while the CEPD model agrees with PIC much better.
Finally, in table I we tabulate the calculated values of ω eff for various angles γ and three different diode scales, i.e. different radii r 0 , α and R, keeping constant ratii R/r 0 = 10 4 and α/r 0 = 0.235. The numerical results confirm the scale invariance of ω shown analytically by eq. (8). We see that ω eff has minute variations -within error margins-for different diode scales. The error margins are obtained from the covariance matrix of the curve fitting to the PIC I − V data.

V. DISCUSSION
The results of figures 4 and 5 demonstrate the basic utility of the CEPD model. If the value of ω eff and the field conversion length χ are available for a certain electrode geometry, the complex problem of calculating the SC suppressed emission current from a 3D emitter is reduced from running full PIC simulations, to evaluating the algebraic formula (4) self-consistently with the emission characteristics (e.g. the Fowler-Nordheim equation).
In order to obtain ω eff theoretically, PIC simulations are required, from which the CEPD model is fitted to the numerical I − V curve. However, this needs to be done only once for a given geometry; then the CEPD model can be used to calculate the emission current for different emission characteristics, such as work function and temperature. This significantly reduces the computational cost of introducing emission routines that include the space-charge effects into more complex simulation techniques (see. e.g. [12]).
Furthermore, ω eff values can be calculated and tabulated for typical emitter geometries becoming readily available for use. The scale invariance of ω eff , demonstrated analytically in section II A and numerically in table I, facilitates this by reducing the free parameters of any geometry. Such a tabulation is out of the scope of this work and shall be given in a forthcoming publication. Alternatively, exactly the same fitting procedure can be performed directly to the experimental I −V curve instead of a numerical one, extracting directly ω eff and giving a direct comparison between experiment and theory.
Finally, a comment is worthy regarding the generality of utilizing a single effective value ω eff for the whole emitter. As mentioned earlier, ω is defined separately for each point on the emitting surface, but if its variance within the area where the bulk of the emission originates is small, the usage of a single value ω eff is sufficient. The same argument holds for the conversion length χ (or the field enhancement factor), allowing for the usage of a single χ value for an emitter, which is typical in the field emission community. This approximation does not hold in case of emitters that have more than one distinct regions contributing significantly to the emission, with different geometrical characteristics each. In this case, the approximation of a single effective value for ω is ex-pected to become invalid, similarly to the invalidation of the standard Fowler-Nordheim theory for such emitters (see. e.g. Ref. [44]).

VI. CONCLUSIONS
We have developed a three-dimensional theoretical model, describing the scaling laws of space charge limited charge emission at high electric fields. Our model generalizes the one-dimensional planar model to be applicable for any geometry, using a geometry-specific correction factor. We validated our model by comparing it to both numerical calculations and existing experimental field emission data, either of which can be used to obtain the geometrical correction factor of our model. We showed that the classical planar model tends to significantly overestimate the space charge effects, whereas our generalized theory is in very good agreement with both numerical calculations and experimental measurements.

ACKNOWLEDGMENTS
The current study was supported by CERN's CLIC K-contract No. 47207461 and the European Union's Horizon 2020 program, under grant No 856705 (ERA Chair "MATTER"). We also acknowledge grants of computer capacity from the Finnish Grid and Cloud Infrastructure (persistent identifier urn:nbn:fi:research-infras-2016072533).
Appendix A: Derivation of the CEPD for the spherical and cylindrical geometries In both the spherical and cylindrical geometries, the solution of the continuity equation can be obtained from the Gauss law, yielding ξ(r) = (R/r) n , where r is the radial coordinate, R is the radius of the emitter, and n = 1, 2 for the cylindrical and spherical case correspondingly. Then the Poisson equation becomes where J s is the current density at the emitter surface, which is uniform due to symmetry. For this problem, it is more convenient to solve the equivalent initial value problem (IVP) rather than the boundary value problem (BVP) addressed previously. Thus, we obtain the potential at radius r as a function of the field on the emitter F . Then this function is inverted in order to obtain F as a function of a fixed potential Φ(r), which corresponds to an applied voltage V = Φ(r) at an electrode residing at r.
For this purpose, we use the reduced variablesr = r/R, φ = Φ/F R and solve the IVP The parameter λ ≡ kJ s √ RF −3/2 , as well as ζ, is indicative of the SC strength. We shall call λ "implicit SC strength", in contrast to the explicit one ζ, because ζ depends on the applied voltage V (or equivalently on the Laplace field F L = V /χ). The latter is the a priori known free variable, unlike F , which is to be obtained.
The solution of (A2) can be expressed in terms of a Volterra integral equation of the second kind as φ(r) = φ 0 (r) + λ r 1 κ(r,r ) where φ 0 (r) is the solution of the Laplace equation (λ = 0), and the integration kernel κ(r,r ) is the Green function for the differential operator in the left hand side of eq. (A2). It is φ 0 = 1 − 1/r, κ = (r −r )/rr (A4) for the spherical case and φ 0 = log(r), κ = log(r/r ) (A5) for the cylindrical one. Eq. (A3) can be solved in the form of an asymptotic power series φ = φ 0 + φ 1 λ + φ 2 λ 2 + · · · (A6) using the Adomian decomposition method. The firstorder term is obtained by inserting φ 0 in the integral and yields φ (S) 1 = (r + 2r log(r))D log(r) −r log(r), for the spherical and cylindrical cases correspondingly. In eq. (A7), D(·) denotes the Dawson function [36]. Now we shall use this to obtain an approximation for the Laplace field F L and the field reduction factor θ = F/F L . Given the potential Φ at distance r, the Laplace field on the emitter is F L = F φ(r)/φ 0 (r). Hence, we can write which gives F L as a function of (F, J, R, r) in the form of a Maclaurin series on J. The inverse function F (F L , J, R, r) can be written in a similar series. Its first order term can be found by evaluating ∂F/∂J, utilizing the implicit function theorem [45]. It yields where ζ ≡ kJ 0 √ Rφ 0 F −3/2 L , as defined in sec. II A. By matching the coefficient of ζ in (A9) and (4), we can obtain the correction factor of the CEPD for the spherical and cylindrical diodes as ω(r) = 3 4 φ 1 (r)φ −3/2 0 (r), yielding eq. (9).