Advanced modeling of field enhanced thermionic emission

Thermionic energy converters (TECs) are class devices, which promise improvements to the efficiency and cost of both small- and large-scale electricity generation.^1,2 A TEC is comprised of a narrowly separated thermionic emitter and an anode. Simple structures are often space-charge limited as operating temperatures produce currents exceeding the Child–Langmuir limit. While in principle these devices are relatively simple to conceive, they face a number of challenges to their practical realization. In particular, achieving the desired device efficiencies requires the development of robust thermionic emitters, able to operate at a range of temperatures, and exhibiting dramatically reduced work functions (on the order of 1–2 eV).³ Novel materials, including complex, layered structures, may be required to meet these requirements.^4–8 

For heterogeneous emitting surfaces, local field enhancement plays a critical role in reducing the effective work function of the emitter.^9,10 In practice, emission characteristics near the space-charge limit are further obscured by the formation of virtual cathodelike modulations.¹¹ Significant efforts have been made to theoretically quantify and measure the field-enhancement factor for different emitter geometries.^12–16 High fidelity simulations are needed to capture the microscale variation in electric fields and current densities along the surface and through the virtual cathode region.

Particle-in-cell codes excel at simulating these kinds of self-consistent, time-dependent particle dynamics but require significant resources to accurately describe the emission process. The Warp framework is a fully 3D particle-in-cell electrostatic code capable of capturing complex time-dependent effects and features advanced capabilities for mesh refinement and subgrid-cell boundaries.¹⁷ Previous work with the Warp code has demonstrated its applicability for modeling intense ion beam sources,¹⁸ in resolving field emission from diamond nanotips,¹⁹ and for simulating thermionic emission in the space-charge limit.²⁰ However, Warp currently employs bulk emission algorithms that cannot be used for arbitrarily shaped thermionic emitters. Moreover, these algorithms use static particle weights, which may result in under- or oversampling of emission surfaces with strongly varying surface fields, as is characteristic of microstructured emitters.

In this paper, we present an improved algorithm for modeling self-consistent, field-enhanced thermionic emission in a finite-difference time-domain (FDTD) scheme using the Warp code. We examine the applicability of subgrid interpolation in conjunction with mesh-refinement for simulating nonconformal, shaped emitters and illustrate convergence properties. We then implement a custom weighting algorithm that permits uniform sampling of emission surfaces with a large variation in the surface electric field. Finally, we demonstrate this algorithm in both the field-enhancement regime and the space-charge-limited regime.

Subgrid boundary applications, or “cut-cell” techniques, have been used in a broad array of uniform-grid codes to improve field accuracy near complex, embedded boundaries.^21,22 Because emission models are extremely sensitive to local electric field variations, we first establish the applicability of subgrids to resolve the expected field enhancement for shaped emitters. We consider different permutations of the subgrid application and the applied stencil used to obtain the field solutions. Table I shows the different subgrid configurations and their corresponding flags in Warp. It should be noted that applying the corrected field calculation when the subgrid is turned off, or when the potential is calculated from the staircase subgrid, is fundamentally flawed. Therefore, these cases are not included in our analysis. For this analysis, we simulate a conical emitter, which can be exactly represented with linear subgrid functions. The emitter has a base radius of $0.2μm$ and a height of $0.2μm$⁠. The gap distance was $1μm$⁠.

We computed the total electrical field on the surface of the cone and the corresponding potential for several different permutations of the settings in Table I. Figure 1 shows the results of these calculations. Because the staircase subgrid artificially enlarges the conductor area, the cut-cell subgrid yields higher fields and potential than the staircase subgrid. The corrected field calculation gives higher and sharper fluctuations in the surface field. These large fluctuations are not desirable when computing the weights of the particles.

We also examined the longitudinal electric field on the axis of the cone for different cut-cell configurations. Figure 2 shows the longitudinal electric field $Ez$ field along the z-axis, which penetrates the center of cone. In the absence of field interpolation and quadratic subgrid interpolation, the location of the emitter tip is not accurately resolved. When the subgrid is turned off but the field interpolation is turned on, the simulation does compute an enhancement factor; however, the field solution is still incorrect. When both the subgrid and field interpolation are used, linear and quadratic interpolations show minor differences.

While the interpolated subgrid calculation provides a more accurate field solution at the tip of the emitter, it produces larger variations in the fields along the emitter surface. Therefore, including the e-field interpolation is not ideal for these shaped emitters. The subgrid potential calculation, however, is clearly necessary for smoothly resolving the fields. During the remainder of our studies, the emitters are simulated with a first order subgrid interpolation for the potential and no subgrid interpolation for the fields.

We next consider the different mesh refinement configurations needed to accurately calculate the fields on the emitter surface and in the simulation domain. The field on the emitter surface is particularly important to ensure accurate calculation of the particle weights in our emission algorithm. While our tool can simulate a range of different emitter tips such as hemispherical, spheroidal, Gaussian, and conical, here we discuss only spheroidal emitters due to the ease of comparison with analytic theory. Figure 3 shows the geometry of the emitter tips used for this study, including locations of the field output and the subgrids for mesh refinement.

First, we examine the convergence of the field solution for different mesh refinement configurations and then compare the field-enhancement factor on the tip of the emitter with analytic theory. Figures 4 and 5 show the field and potential on the surface of a hemispherical emitter, respectively. Good convergence is achieved for a global grid size of $64×64×64$ and the application of two fourfold mesh refinement patches, as indicated by the field error relative to the highest resolution test.

Although the solution converges, a residual nonzero potential remains at the emitter tip surface, indicating the presence of $∼$1% errors that could impact the resulting emission profile. Next, we consider the field along the axis of the emitter. Figures 6 and 7 show the longitudinal electric field and the electrostatic potential along the axis of the simulation domain. We note that maximum mesh refinement is required to resolve the peak electric field at the surface of the emitter tip, located at $z=2.5mm$⁠. While increasing the global mesh density does improve the resolution of the emitter surface, the use of refinement is critical. For surfaces defined by analytical functions, particles are injected exactly on the boundary without knowledge of the simulation mesh. If the emitter is not being accurately resolved, then the field solution used to compute the particle weights will not accurately represent the geometric field enhancement. These results highlight the need for a second refinement patch located at the emitter tip.

Next, we benchmark the field-enhancement factor computed at the tip of the emitter with analytic theory.¹⁵ For this study, we consider a spheroid of varying eccentricity. For a spheroidal emitter, the field-enhancement factor $β$ is given by Eq. (1) in terms of the eccentricity $ξ$⁠. Figure 8 compares the computed field-enhancement factor for a range of spheroidal emitters as a function of the eccentricity compared to Eq. (1). In general, we see very good agreement between Warp and the analytic model. Small errors in the field-enhancement factor are likely due to not totally resolving the tip of the emitter. Because the spherical geometry is defined as a surface of revolution, the emitter tip has a zero transverse derivative, making it impossible to resolve perfectly on a uniform structured mesh,

We have chosen to implement an injection scheme that employs variable particle weights. When simulating field enhanced thermionic emission for complex shapes, this has the advantage of not oversampling the high field region or undersampling the low field region. Weights are computed using the equation

Here, $wi$ is the weight of the $ith$ particle, $Ji$ is the local current density computed by Eq. (3), $Ai$ is the local area where the particle is emitted, $dt$ is the integration step, and $q$ is the fundamental unit of charge,

The use of variable weights also permits the variation of parameters such as the work function (⁠$W(x,y,z)$⁠), temperature (⁠$T(x,y,z)$⁠), material properties (⁠$Ag(x,y,z)$⁠), and local field enhancement (⁠$ΔW(x,y,z)$⁠) along the surface of the emitter. For this work, we are focused primarily on field-enhancement effects; therefore, we assume that the local work function, temperature, and material properties are constant. The local field-enhancement factor is computed in the usual way $ΔW(x,y,z)=q3E(x,y,z)4πε0$⁠, where $E(x,y,z)$ are the fields at the emission site gathered from the solver in Warp.

The local current density is comparatively straightforward to compute, but determining the effective area of the emission site requires care. The simplest method is to sample the surface with a uniform density and apply an effective area based on the number of particles and the surface area of the emitter. For example, a conical emitter would have an effective area given by

This approach is valid given the assumption that local variations in sampling the surface will be averaged out over larger time scales. This method has the advantage of applying a uniform scaling to all particle emission sites based on the area but requires care in determining the emission sites to ensure appropriate sampling. For conical and spherical emitters, this approach is relatively straightforward. However, for more complicated geometries, a generalized approach is needed. We, therefore, divide the surface into a series of domains with a fixed surface area and emit particles from each domain weighted by the local surface area given by

Here, we assume that the local surface is described by analytic functions that are differentiable. For complex surfaces, this integral cannot be computed analytically; however, for a sufficiently small domain size, we can approximate the local surface area using the equation

To validate this emission model in the presence of significant field enhancement, we performed a series of simulations and analyzed the outputs using pseudoempirical methods, as there is no general analytical solution for a shaped emitter in a diode configuration. For the temperature and applied fields in our system, the current through the diode can be modeled by the standard Richardson–Dushamn equation with a spatially dependent field-enhancement factor,

It should be noted that for lower temperatures or higher applied fields, the emission is better described by the general thermal field (GTF) model,^23,24 which has been implemented for shaped emitters in different regimes.^25,26 More exotic structures such as carbon nanotubes may exhibit field emission or hybrid emission properties for applied fields as low as $∼$1 MV/m.^27,28 Accurate representation of such structures using an FDTD particle-in-cell tool is extremely computationally expensive, as resolving the individual tubes requires nanometer-scale computational grids and high order interpolation at the emitting surfaces. For these studies, finite element solvers have been used to resolve the field enhancement and screening factors for individual emitters.²⁹ 

The geometric field-enhancement factor computed from simulation, $β$⁠, should be constant for simulations with the same emission geometry, regardless of the applied field and emitter temperature. Validation of this property provides a convenient check on the consistency of our model. For different simulation parameters, we can compute $β$ by taking the log of the ratio of the current density at the exit of the diode. This approach is described by

Here, we are assuming that the spatial variation in the enhanced field will vary linearly with the applied field and is captured by $β$⁠. Assuming that $β$ does not vary with the temperature or the applied field, the solution is given by

We computed the field-enhancement factor for the spherical emitter operating at four different temperatures for which field enhancement significantly increased the emitted current. Figure 9 shows the field-enhancement factor averaged over five different configurations for the applied voltage. The geometry for this study follows that presented in Fig. 3. For these simulations, a large external magnetic field was applied in order to ensure that the beam is confined transversely. This allows us to calculate $β$ using Eq. (9) and compare with the analytical and numerical results presented in Sec. III. In spite of some numerical noise, the field-enhancement factor remains constant as the temperature and applied field are varied. Additionally, at the tip of the emitter, the field-enhancement factor for each temperature is $∼3$⁠, as expected from the analytical calculations presented earlier. This indicates that our emission model is properly capturing the field-enhancement effects along the emitter tip for different geometries.

Because the field on the spherical tip goes to zero at the edge of the emitter, as seen in Fig. 4, we expect that the field-enhancement factor will also go to zero to account for this variation. Note that for the diode configuration, the applied fields $E1$ and $E0$ are computed as the voltage applied to the anode divided by the gap distance. This variation in the field on the emitter surface will be captured by the $β$ parameter in Eq. (7). The resulting behavior in the output beam is consistent with the field solution on the cathode.

Finally, we consider emission in the regime where self-fields from the particles significantly offset field enhancement and may produce space-charge limited emission. We studied this regime for the hemispherical emitter over a range of temperatures and applied voltages. As with Sec. V, the simulations conducted here included a large longitudinal magnetic field. This confines the beam trajectories to follow a constant transverse position, thereby allowing us to focus only on longitudinal space-charge independently of the strong transverse fields. Figure 10 shows the onset of space-charge limited emission as the temperature is increased for a diode at a constant voltage.

These simulations suggest that we are properly capturing the onset of space-charge limited emission. First, the lack of a current increase from 1600 to 1700 K is a clear indicator of space-charge limited emission. Second, the time-dependent current at the output of the diode is consistent with an emitter that is attempting to push more current than the space-charge limit. This also shows that we are also properly handling time-dependent weights when sampling the fields for the emission of particles. Figure 11 shows the output current as a function of time for a diode operating at 1700 K and several applied voltages.

Here, we see that for higher applied voltages, the height of the peak transmitted current is slightly lower. This is consistent with space-charge limited transport for a diode operating at different applied voltages. Furthermore, as with Sec. V, we can examine the geometric field-enhancement factor in the space-charge limited configuration. For a constant geometry, we expect the output current to scale as $J=g(βe)2J0$⁠.³⁰ Here, $g(βe)2$ is a geometric factor and $J0$ is the space-charge limited current for the diode configuration without field enhancement. We computed the geometric factor $g(βe)2$ for a range of voltages and a cathode temperature of 1700 K (Fig. 12).

Here, we see consistent solutions for the geometric field-enhancement factor in the space-charge limited regime, further indicating that our models are capturing these behaviors. Small fluctuations between different parameters can be attributed to noise in the field solution and in the particle statistics.

In this paper, we have demonstrated an improved algorithm for modeling self-consistent, field-enhanced thermionic emission from complex surfaces in a finite-difference time-domain scheme using the Warp code. We validated solutions to the electrostatic potential and the corresponding electric field using a linear cut-cell representation of the emission boundary and illustrated the use of mesh refinement for better resolving surface fields from complex geometries. We then implemented a custom emission model that employs variable particle weights to uniformly sample emission surfaces with a large range of surface fields. We validated this emission model by numerically computing the geometric field-enhancement factor for a range of configurations. Finally, we verified that space-charge effects are properly captured when cathode temperatures are large relative to the applied voltage.

Our results lay the foundation for the use of Warp to simulate and optimize complex vacuum electronic devices. The use of the described emission models enables the simulation of a wide range of devices with variable geometries, work functions, and temperatures. These improved capabilities greatly expand the range of devices that can be modeled in Warp. Bringing improved high fidelity modeling and optimization to the broader vacuum electronics community has the potential for a profound impact on how devices are studied and designed in the future.

This material is based upon work supported by the U.S. Department of Energy (DOE), Office of Science, Office of Advanced Scientific Computing Research under Award No. DE-SC0017162.