title: General scaling laws of space charge effects in field emission
authors: Andreas Kyritsakis Mihkel Veske Flyura Djurabekova
text: 
  
  Sample title]General scaling laws of space charge effects in field emission
  
  
  akyritsos1@gmail.com
  Helsinki Institute of Physics and Deparment of Physics, University of Helsinki, PO Box 43 (Pietari Kalmin katu 2), 00014 Helsinki, Finland
  Institute of Technology, University of Tartu, Nooruse 1, 50411 Tartu, Estonia
  
  
  Helsinki Institute of Physics and Deparment of Physics, University of Helsinki, PO Box 43 (Pietari Kalmin katu 2), 00014 Helsinki, Finland
  
  
  
  
  
  Helsinki Institute of Physics and Deparment of Physics, University of Helsinki, PO Box 43 (Pietari Kalmin katu 2), 00014 Helsinki, Finland
  National Research Nuclear University MEPhI, Kashirskoye sh. 31, 115409 Moscow, Russia
  
  
  
  
  
  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.
  
  
  
  
  
  [
      F. Djurabekova
      December 30, 2023
  =====================
  
  
  
  
  
  § 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 <cit.>, 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 <cit.> and is of paramount importance for the understanding of the ignition of vacuum arcs <cit.>.
  
  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) <cit.>, molecular dynamics <cit.>, or the point charge method <cit.>.
  
  Analytical solutions of the SC problem exist only for 1D geometries, where the continuity equation has a simple solution.
  For instance, Child <cit.> and Langmuir <cit.> 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 <cit.> and spherical <cit.> electrodes. 
  In a recent work <cit.>, 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. <cit.>.
  This solution has been widely used as a reference model to estimate the SC effects on field emission <cit.>, 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. <cit.>, i.e. by assigning the real values of voltage and cathode field to a planar diode.
  
  However, PIC simulations <cit.> 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 <ref> 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 geometry-dependent correction factor ω.
  In section <ref> we derive ω for the spherical and cylindrical diode problems.
  In section <ref> 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 <ref>, showing that the correction factor can be obtained for a certain geometry by fitting to either numerical calculations or experimental data.
  
  
  
  § THEORY
  
  
  
  
   §.§ General formulation
  
  
  The standard formulation of the space charge problem adopted already since Child's work <cit.>, assumes a continuous charge density distribution of the emitted charge ρ(𝐫) and zero initial velocity for the emitted charged particles.
  Some studies have considered the case of non-zero initial kinetic energy <cit.>, 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
  
      ∇^2 Φ =kJ(𝐫)Φ^-1/2,
  
  with boundary conditions Φ = 0 at the emitter and Φ = V at the collector electrodes.
  Here Φ is the electrostatic potential, k=ϵ_0^-1√(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(𝐫) obeys the continuity equation ∇·𝐉 = 0, with boundary condition 𝐉(𝐫_s) = J_s(𝐫_s)n̂(𝐫_s), where J_s(𝐫_s) is the emitted current density and n̂(𝐫_s) is the normal unit vector at the emitting surface point 𝐫_s. 
  In order to obtain 𝐉 and Φ, the above equations have to be solved self-consistently, along with the surface emission laws that give J_s(𝐫_s) as a function of the local electric field F(𝐫_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 <cit.>
  
      6(kJ)^2d-F^3 = (2kJ√(V)-F^2 ) √(4kJ√(V)+F^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. (<ref>) needs to be solved self-consistently with the law that gives J(F), e.g. the Fowler-Nordheim (FN) law <cit.> or its modern generalizations <cit.>.
  Barbour et. al. <cit.> 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 <cit.>.
  
  Eq. (<ref>) can be simplified by introducing the reduced dimensionless variables <cit.> θ≡ Fd/V=F/F_L and ζ≡ kJd^2V^-3/2 = kJd^1/2F_L^-3/2 = kJV^1/2 / F_L^2, yielding
  
      3θ^2(1- θ) = ζ(4 - 9ζ) .
  
  In (<ref>), θ 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 (<ref>) that θ reduces from 1 to 0 (F from F_L to 0) as ζ increases from 0 to 4/9, where the Child law <cit.> limit occurs.
  Eq. (<ref>) can be solved analytically <cit.>; yet, it is more convenient to express the physical solution in a perturbation series around ζ=0, yielding
  
      θ = 1- 4/3ζ - 5/9ζ^2 - 16/27ζ^4 + ⋯.
    
  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^12A/m^2) <cit.>.
  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 𝐫_s at the emission surface, with local field F.
  We shall express the Poisson equation in terms of the reduced variables ϕ≡Φ / V, J(𝐫)=J_s ξ(𝐫), where J_s = J(𝐫_s) and ξ(𝐫) is a unitless variable expressing the variation of the current density from 𝐫_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 ξ(𝐫) depends only on the emitter geometry and not on J_s.
  Finally, we use the reduced space coordinate 𝐫̃ = 𝐫 / χ, where χ is the “conversion length” χ≡ V / F_L (F_L is the Laplace (J=0) field at 𝐫_s).
  Equation (<ref>) becomes
  
      ∇̃^2 ϕ(𝐫̃) = ( kJ_sχ^2/V^3/2) ξ(𝐫̃)/√(ϕ(𝐫̃)),
  
  where ∇̃ denotes derivatives with respect to 𝐫̃. 
  Note that the reduced position r̃ is scale invariant, i.e. it remains unchanged under a geometrical scaling r→ a r, 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 χ^2V^-3/2=kJ_s V^1/2/F_L^2 with the space charge strength of the planar diode.
  The classical planar diode equivalence <cit.> corresponds to inserting the above ζ to eq. (<ref>), or equivalently substituting d with χ in eq. (<ref>).
  Now we shall derive an asymptotic expansion similar to eq. (<ref>) 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. (<ref>) in its integral form, utilizing the Green's function for Dirichlet problems <cit.> G(𝐫̃, 𝐫̃'), i.e. the solution of the boundary value problem (BVP) ∇̃_𝐫^2 G(𝐫̃, 𝐫̃') = δ(𝐫̃- 𝐫̃')  on Ω,  G(𝐫̃, 𝐫̃') = 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
  
      ϕ(𝐫̃) = ϕ_0(𝐫̃) + ζ∫_ΩG(𝐫̃, 𝐫̃') ξ(𝐫̃')/√(ϕ(𝐫̃')) d^3𝐫̃' ,
  
  where ϕ_0(𝐫̃) is the Laplace solution (ζ=0).
  For ζ≪ 1, ϕ(𝐫̃) can be expanded in an asymptotic series using the Adomian decomposition method <cit.>, as
  
      ϕ(𝐫̃) = ϕ_0(𝐫̃) + ζ∫_ΩG(𝐫̃, 𝐫̃') ξ(𝐫̃')/√(ϕ_0(𝐫̃')) d^3𝐫̃' + O(ζ^2) .
  
  By taking the gradient we obtain a similar expansion for the field reduction factor θ = F / F_L
  
      θ   = 1 - 4/3ωζ + O(ζ^2) , with
      ω   ≡3/4∫_Ω.∇̃_𝐫 G(𝐫̃, 𝐫̃')|_𝐫̃_sξ(𝐫̃')/√(ϕ_0(𝐫̃')) d^3𝐫̃' .
  
  
  The above equation demonstrates the inadequacy of the standard planar diode equivalence of Refs <cit.>.
  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 𝐫_s, we define the corrected equivalent planar diode model, for which the SC strength is ζ_c = ω kJ_sχ^2V^-3/2=ω kJ_s√(χ)/F_L^3/2.
  Practically, the CEPD model can be used by replacing ζ = ζ_c in eq. (<ref>), or using the CEPD gap distance d_c = ω^2 χ and voltage V_c = ω^2V in eq. (<ref>) and solving for F.
  Alternatively, one can substitute directly V and d=χ in  (<ref>), as in Refs <cit.>, 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. (<ref>) 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 <cit.>.
  Here we define and calculate it rigorously.
  ω can be formally calculated for a given emitter geometry and point 𝐫_s using (<ref>), 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.
  
  
  
   §.§ 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 <cit.>. 
  Langmuir and Blodgett <cit.> solved the specific SC problem for both geometries.
  Furthermore, for the spherical case, which has been recently used to model sharp field emitters <cit.>,  Aizenberg <cit.> 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
  
      ω^(S)    = 3/4(2- 1/r̃) log( √(r̃) + √(r̃-1)) - √(1- 1/r̃)/( 1 - 1/r̃)^3/2,
      ω^(C)   =  3/4(r̃+2 r̃log (r̃)) D(√(log (r̃)))-r̃√(log (r̃))/[log(r̃) ]^3/2,
  
  for the spherical and cylindrical diodes correspondingly, where r̃ is the ratio between the emitter and the collector radii and D(·) denotes the Dawson function <cit.>.
  The derivation is given in appendix <ref>.
  
  It is worth noting that for r̃→ 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, for r̃→∞, 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 increasing r̃.
  
  
  
  
  
  In order to validate the above results, we solved the Poisson equation numerically by the Runge-Kutta method <cit.>.
  The results are shown in figure <ref> where we plot the field reduction factor θ vs both the SC strength ζ (dashed lines) and the corrected one ζ_c = ωζ (solid lines), for various values of r̃, 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 ratio r̃.
  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. (<ref>) (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 with r̃, with the increase being much faster for the cylindrical diode.
  
  
  
  § METHODS
   
  
  
   §.§ 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 <cit.>, which has been recently enhanced with PIC capabilities <cit.>.
  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 sub-picosecond timescale <cit.>.
  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 <cit.>.
  
  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 <ref>).
  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 <cit.>, which evaluates the appropriate emission formulas.
  Then charged superparticles (SPs) are injected from the surface into the vacuum.
  Unlike our previous work <cit.>, 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. <cit.>, we use the explicit leapfrog method integration scheme <cit.>.
  However, since here we are integrating only the electron paths without a concurrent 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 <ref>, 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. <cit.> 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. (<ref>) 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.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ 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. <cit.>, an electrochemically etched tungsten cathode was utilized to take I-V measurements at both low and high currents. 
  The cathode 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 <cit.>, including the one of Ref. <cit.>, are generally described by the general sphere-on-a-cone (SOC) model, developed by the same group <cit.>. 
  
  
  
  
  
  
  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+1r^-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 <cit.> 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. (<ref>).
  
  Figure <ref> 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. <cit.> 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. <cit.> 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. <cit.>) β≡χ^-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. <cit.>.
  
  
  
  
  § RESULTS
   
  
  Given the geometrical model described in section <ref>, we used the computational method described in section <ref> 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 <ref> 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 <ref>, varies on the emitting surface, mainly due to the variation of the current density distribution ξ(𝐫̃). 
  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 <ref>, we compare the experimental data of Ref. <cit.> (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. <cit.>.
  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. <cit.>. 
  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 effective emission area due to the increased non-uniformity of the work function, evident from the corresponding micrograph of figure 5 in Ref. <cit.>.
  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 <ref>, 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 <cit.>, 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. 
  
  
  However, this is rather a coincidence, specific to this particular geometry.
  In figure <ref> 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 <ref>. 
  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 <ref> 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. (<ref>). 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.
  
  
  
  
  
  
  § DISCUSSION
  
  The results of figures <ref> and <ref> 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 (<ref>) 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. <cit.>). 
  
  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 <ref> and numerically in table <ref>, 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 expected to become invalid, similarly to the invalidation of the standard Fowler-Nordheim theory for such emitters (see. e.g. Ref. <cit.>).
  
  
  
  
  § 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).
  
  
  
  
  
  
  
  § 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
  
      d/dr( r^n dΦ/dr) = kJ_sR^n/√(Φ),
  
  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 variables r̃ = r/R, ϕ = Φ/FR and solve the IVP
  
      d/dr̃( r̃^n dϕ/dr̃)     = kJ_s√(R)/F^3/21/√(ϕ) = λ/√(ϕ),
      . d ϕ/dr̃|_r̃=1    = 1 .
  
  The parameter λ≡ kJ_s√(R)F^-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 (<ref>) can be expressed in terms of a Volterra integral equation of the second kind as
  
      ϕ(r̃) = ϕ_0(r̃) + λ∫_1^r̃κ(r̃,r̃')/√(ϕ(r̃'))dr̃' ,
  
  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. (<ref>). 
  It is
  
      ϕ_0 = 1 - 1/r̃, κ = (r̃ - r̃')/ r̃r̃'
  
  for the spherical case and 
  
      ϕ_0 = log(r̃), κ = log(r̃/r̃')
  
  for the cylindrical one.
  Eq. (<ref>) can be solved in the form of an asymptotic power series
  
      ϕ = ϕ_0 + ϕ_1 λ + ϕ_2 λ^2 + ⋯
  
  using the Adomian decomposition method.
  The first-order term is obtained by inserting ϕ_0 in the integral and yields
  
      ϕ_1^(S)    = (2- 1/r̃) log( √(r̃) + √(r̃-1)) - √(1- 1/r̃),
      ϕ_1^(C)   = (r̃+2 r log (r̃)) D(√(log (r̃)))-r̃√(log (r̃)),
  
  for the spherical and cylindrical cases correspondingly.
  In eq. (<ref>), D(·) denotes the Dawson function <cit.>.
  
  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
  
      F_L= F (1 + λϕ_1/ϕ_0 + O(λ^2) ),
  
  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 <cit.>.
  It yields 
  
      F = F_L ( 1 - ζϕ_1/ϕ_0^3/2 + O(ζ^2) ) ,
  
  where ζ≡ kJ_0√(Rϕ_0)F_L^-3/2, as defined in sec. <ref>. By matching the coefficient of ζ in (<ref>) and (<ref>), we can obtain the correction factor of the CEPD for the spherical and cylindrical diodes as ω(r̃) = 3/4ϕ_1(r̃)ϕ_0^-3/2(r̃), yielding eq. (<ref>).
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   