Thermionic Emission in a Planar Diode

In this benchmark model, electrons are released via thermionic emission from a hot cathode in a plane parallel vacuum diode. Because the thermal electrons contribute to the space charge density in the region between the cathode and the anode, a potential barrier forms close to the cathode. This potential barrier retards the motion of electrons toward the anode so that electrons of lower energy are repelled back toward the cathode and only electrons with higher energy are transmitted. As a result, the transmitted current in the diode is only a fraction of the total thermionic current.

The Charged Particle Tracing interface is used to model thermionic emission from a cathode with a specified temperature and work function. A dedicated study for modeling bidirectionally coupled particle field interactions is used to include the contribution of the thermal electrons when computing the electric potential in the diode. The transmitted current and the potential in the diode are compared to the analytic solution of the Langmuir-Fry model.

The phenomenon of thermionic emission can be subdivided into different regimes based on the temperature and potential difference. In both regimes, a population of electrons following the classical Maxwell-Boltzmann energy distribution are emitted from a cathode at a fixed temperature. The emitted electrons have a velocity direction distribution following Lambert’s cosine law. For a fixed temperature, as the electric potential difference between the cathode and anode is increased, the collected current at the anode increases and will eventually saturate. This special condition is called the thermally limited regime; the current at the anode depends only on the work function and temperature of the cathode and is unaffected by further increases in the potential difference. Before reaching this regime, at lower potential difference, the space charge in the region adjacent to the cathode creates a region of negative potential. This potential barrier repels some of the thermally emitted electrons back to the cathode surface, so the current at the anode is less than the emitted current at the cathode or saturation current. In this case, the electron emission is said to be in the space charge limited regime.

Space charge limited emission of thermal electrons in a plane parallel vacuum diode is one of the only problems in nonlaminar charged particle beam physics that has an exact analytic solution, making it a useful benchmarking tool. This solution is often called the Langmuir-Fry model, referring to the authors of two frequently cited early publications in this area (1,2).

Model Definition

The planar diode consists of two semi-infinite parallel electrodes separated by a distance d = 1 mm. The model geometry is three-dimensional but can be understood as quasi-1D because the potential and current are uniform in the x- and z-directions. The geometry of the diode and the functional form of the electric potential distribution in the space charge limited regime are illustrated in Figure 1. The cathode is grounded (Vc = 0) and the anode is maintained at a fixed potential of Va = 100 V. The cathode is maintained at a uniform temperature of Tc = 2500 K and its work function is Φm = 4.5 V.

Note that Figure 1 is not necessarily drawn to scale. The negative potential barrier is often extremely narrow compared to the diode width d.

Figure 1: The planar diode in the space charge limited regime. The black solid curve shows the electric potential between the electrodes, which reaches a local minimum at the potential barrier before increasing closer to the anode.

The geometry is oriented so that the normal vector between the two parallel electrodes points in the positive y-direction. From the symmetry of the semi-infinite planar diode, the transverse electric field components are zero (Ex = Ez = 0), and the components of the electron velocity in these directions remain unchanged after they leave the cathode (vx = vx0, vz = vz0). For this reason, a mesh is created with only a single element in the transverse directions.

This model uses the dedicated Particle Field Interaction, Non-Relativistic interface, which creates the following interfaces and couplings:

•

An instance of the Electrostatics interface to compute the electric potential,

•

An instance of the Charged Particle Tracing interface to exert forces on the electrons and compute their trajectories, and

•

An Multiphysics coupling to compute the effect of the particles on the space charge density distribution within the modeling domain.

Although the electrons always exist at a number of discrete positions, the charge density must be represented as a field with finite values. The node defines a contribution from each model particle to the mesh element that the particle occupies, allowing their charge to be converted to a spatial density.

A further advantage of using a mesh with only one element in each cross-section parallel to the xz-plane is that there is no possibility of non-physical transverse electric field components arising from discretization of the charge density and discrete sampling of the electrons in velocity space.

The Thermionic Emission Feature

In this example, electrons are released from the cathode using the dedicated release feature. This feature pseudorandomly samples the initial velocity components of the emitted electrons so that they follow a Maxwell-Boltzmann energy distribution and a Lambertian distribution of release angles with respect to the surface normal.

The initial particle velocity components are

where vn is the velocity component normal to the cathode surface and vt1 and vt2 are the two orthogonal velocity components parallel to the surface. The azimuthal angle φ (SI unit: rad) and the polar angle θ (SI unit: rad) follow the probability distribution function

That is, φ has a uniform probability in the interval [0, 2π] and θ has a Lambertian probability distribution in the interval [0, π/2].

The velocity magnitude follows a probability distribution function

where W (dimensionless) is the normalized particle kinetic energy, defined as

where

•

me = 9.10938188 × 10-31 kg is the electron mass,

•

kB = 1.3806488 × 10-23 J/K is the Boltzmann constant, and

•

Tc (SI unit: K) is the cathode temperature.

For details on how the initial particle velocity components are sampled from these distributions, see the Particle Tracing Module User’s Guide.

Bidirectionally Coupled Particles and Fields

Just as the electrons contribute to the space charge density which affects the electric potential, changes in the potential modify the electron trajectories. In this way the particle trajectories and field are bidirectionally coupled. Computing a self-consistent solution, in which the effects of each part of the solution on the other are taken into account, requires the following steps:

Compute the particle trajectories, assuming no space charge effects are present, using a . From these trajectories, compute the space charge density using the node. Because there is no potential barrier during the first iteration, all of the particles should reach the anode and the transmitted current should equal the thermionic current at the cathode.

Compute the electric potential due to the space charge density of electrons in the diode, using a . The boundary conditions for the Electrostatics interface are the specified potentials at the cathode and the anode.

Use the electric potential calculated in step 2 to compute the perturbed particle trajectories. Recalculate the space charge density using these perturbed trajectories.

Repeat steps 2 and 3 until either a specified number of iterations has been reached or until the change in the transmitted thermionic current between successive iterations is less than a specified tolerance.

This iterative procedure is automatically set up when using the study step, a specialized study step intended specifically for modeling bidirectional couplings between time-dependent particle positions and stationary fields.

The Langmuir-Fry Model

The following analytic solution for the electric potential in a plane parallel vacuum diode in the space charge limited regime is based on Ref. 3.

The maximum total current density Jth (SI unit: A/m2) that can be drawn from a cathode has been predicted theoretically and shown experimentally to follow the Richardson equation,

(1)

where

•

A∗ (SI unit: A/(m2K2)) is the effective Richardson constant,

•

Tc (SI unit: K) is the cathode temperature,

•

e = 1.602176565 × 10-19 C is the fundamental charge, and

•

kB = 1.3806488 × 10-23 J/K is the Boltzmann constant, and

•

Φ (SI unit: V) is the cathode work function.

Jth is also called the saturation current. If the voltage saturation has not been reached, and thus the current at the anode is less than Jth, then there must be a potential minimum between the cathode and anode that repels electrons of low energy back toward the cathode as illustrated in Figure 1.

It is convenient to express the charge density and current density in terms of the number density of particles in 6-dimensional phase space n(x,v) (SI unit: s3/m6), where x (SI unit: m) is the position vector with components x, y, and z; and v (SI unit: m/s) is the velocity vector with components vx, vy, and vz.

The electric potential V (SI unit: V) in the vacuum diode follows Poisson’s equation,

(2)

where ε0 = 8.854187817 × 10-12 F/m is the permittivity of vacuum and ρ (SI unit: C/m3) is the charge density of the electrons. In terms of the phase space number density, the charge density is

(3)

In terms of the phase space number density, the current density of electrons emitted at the cathode is

(4)

Note that the range of integration starts from 0 in the y-direction since the velocity of the electrons escaping the cathode can only be positive in this direction.

Also note that the sign conventions have been chosen here so that Jth is positive for the electron emission, even though the net flow of charge is in the -y-direction. This is consistent with the sign convention used in Equation 1 for the saturation current.

The phase space number density is assumed to follow a Gaussian distribution in velocity space,

(5)

where the coefficient F (SI unit: s/m) is given by

and where me = 9.10938188 × 10-31 kg is the electron mass.

The constant C (SI unit: s3/m6) can be expressed in terms of the saturation current by substituting Equation 5 into Equation 4 and integrating, which yields

or alternatively

Thus, in terms of the saturation current the phase space number density is

(6)

Liouville’s theorem states that the phase space density is invariant along a particle’s trajectory,

where the subscripts 0 and 1 indicate the particle’s coordinates in phase space at two different times. Due to the symmetry of the planar diode, the phase space density is uniform in the x- and z-directions, so these arguments can be dropped. Furthermore, since the electric potential in the semi-infinite planar diode only varies in the y-direction, the transverse velocity components are constant along the particle’s phase space trajectory. Thus, Liouville’s theorem can be simplified to

Taking x0 and v0 to be the phase space coordinates at the cathode where y = 0, this can be further simplified to

Substituting this result into Equation 5 yields an expression for the phase space number density at the cathode in terms of the initial velocity components at the cathode v0x, v0y, and v0z,

In the space charge limited regime, electrons are either transmitted to the anode or repelled back to the cathode depending on whether they have sufficiently high velocity in the y-direction to bypass the negative potential barrier close to the cathode. From conservation of energy for nonrelativistic electrons, the y-component of the velocity at any point along an electron’s trajectory is a function of the electric potential V at that point and the electron’s initial velocity in the y-direction,

(7)

where η is the charge-to-mass ratio of the electrons,

The electric potential is measured with respect to the cathode; that is, Vc = 0. The electric potential reaches a minimum value at y = ym (0 < ym < d, where d (SI unit: m) is the distance from the cathode to the anode) where the following conditions apply:

(8)

Because V = 0 at the cathode and the electrons are negatively charged, Vm must be negative. The normal component of the electric field Ey is positive between the cathode and the semi-infinite plane located at y = ym (region I in Figure 1). This positive electric field in region I decelerates the electrons. For y > ym (region II in Figure 1), Ey is negative and the electrons are accelerated toward the anode. In order to enter region II, the emitted electrons must have sufficiently high velocity in the y-direction to bypass the potential barrier at ym. Substituting vy = 0 and V = Vm into Equation 7 yields the minimum initial electron velocity needed to reach the location of the minimum potential,

(9)

In region I, electrons can move in both the positive and negative y-directions. An electron will cross a plane at y (0 < y < ym) with electric potential V (0 > V > Vm) exactly once (going in the positive y-direction only) if and only if

since the electron then has sufficient energy to bypass the barrier and doesn’t return to the cathode.

An electron will cross the plane exactly twice, first going in the positive y-direction and then returning to the cathode in the negative y-direction, if and only if

since in this case the particle has sufficient energy to bypass the plane but insufficient energy to bypass the entire potential barrier, so it eventually gets repelled back to the cathode.

In region II, a particle passes a plane at y (y > ym) once if and only if

since, once a particle has bypassed the potential barrier, it will definitely reach the anode.

The minimum possible value of the y-component of the particle velocity vmy at any location is therefore

(10)

This result will be used in subsequent integrations over velocity space.

Substituting Equation 3 for the space charge density into Equation 2 for the electric potential, taking into account the invariance of the electric potential in the transverse directions, and adjusting the lower limit of integration based on Equation 10 yields the following 1D Poisson equation for the electric potential V at any position y in the diode:

(11)

This integro-differential equation is subject to the boundary conditions V(0) = 0 and V(d) = Va. The lower limit of integration on y is given by Equation 10.

This integro-differential equation can be solved to compute the location and magnitude of the minimum potential. It is convenient to replace the electric potential with the dimensionless variable Y,

(12)

Algebraic manipulation of Equation 11 (see for example Ref. 2) leads to

(13)

where Yc is the value of the nondimensionalized potential Y at the cathode,

Ya is the value of Y at the anode,

or in terms of Yc,

and the dimensionless functions QI(Y) and QII(Y) have been defined as follows:

Thus QI(Y) is used when integrating over values of the nondimensionalized potential in region I, and similarly QII(Y) is used when integrating over values of the nondimensionalized potential in region II.

The current density at the anode is

(14)

where instead of being zero as in Equation 4, the lower limit of integration v0my is the minimum velocity needed for the emitted particles to cross the potential barrier and thus reach the anode. Substituting the number density from Equation 6 yields

(15)

This integral can be rearranged as follows:

Integrating the above expression, and noting that, by definition,

yields a simplified expression for the anode current in terms of the known thermionic current and the nondimensionalized potential at the cathode,

This leads to definitions of Yc and Ya in terms of the cathode and anode currents,

(16)

(17)

Thus, given the electron properties, thermionic current, and potential difference between the cathode and the anode, it is possible to solve numerically for the anode current Ja by substituting Equation 16 and Equation 17 into Equation 13 and then evaluating the integrals over Y.

Any value of Y corresponds to a unique value of the electric potential greater than or equal to the minimum potential Vm. This value of the potential may correspond to a value of y on one or both sides of the potential minimum. The corresponding value of y in region I, denoted yI, is

(18)

The corresponding value of y in region II, denoted yII, is

(19)

Similarly it is possible to compute the potential at a specified y-coordinate. One viable approach is to integrate Equation 18 and Equation 19 for a list of Y-values, then generate a pair of interpolation functions giving yI and yII as approximate functions of Y. These interpolation functions can be inverted and combined to display the electric potential as a function of the y-coordinate in the diode.

Results and Discussion

Figure 2 shows the electric potential profile near the cathode for both the particle tracing simulation and the analytical solution. The curve was created by defining an interpolation function containing numerical integrations of the results of the previous section.

Figure 2: Electric potential next to the cathode. The black solid curve represents the electric potential calculated using the numerical integration. The solid blue curve shows the electric potential obtained from the thermionic emission feature.

Table 1 shows a comparison of the anode current density Ja, minimum potential Vm and minimum potential location ym obtained from the model and from the analytical solution.

The values in the column may change slightly when the model is run on different architectures because the initial velocity of each electron is generated pseudorandomly.

Table 1: CompariSon between the model and the analytical solution.
Parameter	Model	Reference	Relative error (%)
Ja (A/cm2)	-0.26932	-0.25356	6.2150
Vm (mV)	-179.97	-179.37	0.33268
ym (μm)	11.788	11.412	3.2946

References

1. T.C. Fry, “The thermionic current between parallel plane electrodes; velocities of emission distributed according to Maxwell's law,” Physical Review, vol. 17, no. 4, p. 441, 1921.

2. I. Langmuir, “The effect of space charge and initial velocities on the potential distribution and thermionic current between parallel plane electrodes,” Physical Review, vol. 21, no. 4, 419, 1923.

3. P. T. Kirstein, G. S. Kino, and W. E. Waters, “Space charge flow,” in Physical and Quantum Electronics Series, McGraw-Hill, pp. 265–275, 1967.

Particle_Tracing_Module/Charged_Particle_Tracing/planar_diode

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Geometry 1

Load the parameters for the geometry and physics setup from a file. This file also includes the expected values of the anode current, minimum electric potential, and location of the minimum potential, which are also given in the column in Table 1.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file planar_diode_parameters.txt.

Geometry 1

Block 1 (blk1)

In the toolbar, click

In the window for , locate the section.

In the text field, type rc.

In the text field, type d.

In the text field, type rc.

Locate the section. In the text field, type -rc/2.

In the text field, type -rc/2.

Add a thin layer adjacent to the cathode. This will be used to plot the electric potential in a small region surrounding the potential barrier, making the results easier to visualize.

Click to expand the section. Find the subsection. Clear the check box.

Select the check box.

In the table, enter the following settings:


Layer name	Thickness (m)
Layer 1	d/20

The section The Langmuir-Fry Model includes a derivation of the electric potential in the planar diode that can be integrated numerically. Load this numeric data from the reference file and use it to define an interpolation function.

Definitions

Interpolation 1 (int1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Click .

Browse to the model’s Application Libraries folder and double-click the file planar_diode_reference.txt.

Click .

In the text field, type reference.

Locate the section. In the text field, type um.

In the text field, type mV.

Define a nonlocal minimum coupling at one of the edges in the thin layer adjacent to the cathode. This coupling will be used to compute the value and location of the minimum electric potential.

Minimum 1 (minop1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Select Edge 4 only.

Materials

Material 1 (mat1)

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	1	1	Basic

Electrostatics (es)

Set up the boundary conditions for the calculation of the electric potential.

Ground 1

In the window, under right-click and choose .

Select Boundary 2 only.

Electric Potential 1

In the toolbar, click

and choose .

Select Boundary 9 only.

In the window for , locate the section.

In the V0 text field, type Va.

Charged Particle Tracing (cpt)

In the window, under click .

Thermionic Emission 1

In the toolbar, click

and choose .

Select Boundary 2 only.

In the window for , locate the section.

In the T text field, type Tc.

In the Φ text field, type Phim.

In the A* text field, type A0.

In the N text field, type N.

Locate the section. In the Nvel text field, type Nvel.

From the list, choose . Sampling at uniform speed intervals causes a disproportionately large number of model particles to have above-average speed. However, since the faster particles are more likely to cross the barrier and thus have more significant contributions to the charge density near the anode, it is reasonable to devote more degrees of freedom to them. This disproportionate sampling from the probability distribution function is offset by reducing the contributions of these energetic particles to the space charge density in the domain.

In the n text field, type 20.

Particle Properties 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Apply the boundary condition to the sides of the modeling domain. The condition acts like a specularly reflecting boundary. The physical interpretation of this is that, on average, for every electron that leaves through the sides of the modeling domain, another electron enters the domain at the same time. This is consistent with the interpretation of the domain as a finite region of a semi-infinite diode.

Symmetry 1

In the toolbar, click

and choose .

Select Boundaries 1, 3–5, 7, 8, 10, and 11 only.

Use the feature to create a global variable for the total current at the anode.

Particle Counter 1

In the toolbar, click

and choose .

Select Boundary 9 only.

Electric Force 1

In the window, click .

In the window for , locate the section.

From the E list, choose .

Locate the section. Select the check box.

To make the solution more robust and reproducible, modify the node so that the space charge density used to compute the electric potential is the cumulative average over successive iterations of the solver sequence.

Multiphysics

Electric Particle Field Interaction 1 (epfi1)

In the window, under click .

In the window for , locate the section.

Select the check box.

In the β text field, type 0. A value of zero means that the space charge is not ramped up for the purpose of computing the electric potential. Instead the full space charge contribution will be used even in the first iteration.

From the list, choose . Weighting the iterations in an arithmetic sequence minimizes the impact of the earlier iterations, which generally use only crude initial guesses for the potential distribution and will have larger relative error compared to later iterations.

The planar diode is essentially a one-dimensional model. Mesh the transverse directions of the geometry with only one mesh element to reduce the number of degrees of freedom. An additional benefit is that no spurious transverse electric forces can arise due to the discretization of the particle distribution in velocity space.

Mesh 1

Mapped 1

In the toolbar, click

and choose .

Select Boundary 2 only.

Distribution 1

Right-click and choose .

In the window for , locate the section.

In the text field, type 1.

Select Edges 1, 3, 5, and 14 only.

Swept 1

In the toolbar, click

Distribution 1

Right-click and choose .

In the window for , locate the section.

Click

Select Domain 1 only.

Locate the section. From the list, choose .

In the text field, type 80.

In the text field, type 10.

Select the check box.

Distribution 2

In the window, right-click and choose .

In the window for , locate the section.

Click

Select Domain 2 only.

Locate the section. From the list, choose .

In the text field, type 100.

In the text field, type 10.

Click

. The resulting swept mesh should be extremely fine near the cathode and coarser at greater distances from the cathode.

Increase the number of iterations in the solver sequence. As more iterations are taken, it becomes easier to obtain a more self-consistent solution.

Study 1

Step 1: Bidirectionally Coupled Particle Tracing

In the window, under click .

In the window for , locate the section.

From the list, choose .

In the text field, type 0 0.1.

Locate the section. In the text field, type 50.

In the toolbar, click

Results

Use the dataset to plot the electric potential in a narrow region adjacent to the cathode.

Edge 3D 1

In the toolbar, click

and choose .

Select Edge 4 only.

Electric Potential Minimum

In the toolbar, click

In the window for , type Electric Potential Minimum in the text field.

Locate the section. From the list, choose .

From the list, choose .

Locate the section. Select the check box.

In the associated text field, type Electric potential (mV).

Locate the section. From the list, choose .

Line Graph 1

Right-click and choose .

Use the default expression, which is V, the electric potential.

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the text field, type y.

From the list, choose .

Click to expand the section. From the list, choose .

Select the check box. In the table, enter Particle tracing.

Line Graph 2

In the window, right-click and choose .

In the window for , locate the section.

In the text field, type reference(y).

From the list, choose .

Locate the section. From the list, choose .

In the text field, type y.

From the list, choose .

Click to expand the section. From the list, choose .

Locate the section. Select the check box.

From the list, choose . In the table, enter Reference.

In the toolbar, click

. Compare the resulting plot to Figure 2.

Global Evaluation 1

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Locate the section. In the table, enter the following settings:


Expression	Unit	Description
cpt.pcnt1.It	A	Transmitted current, Particle Counter 1
minop1(V,y)	um	Location of the potential minimum
minop1(V)	mV	Potential minimum

Compare the resulting values to the column in Table 1. The values may differ slightly due to the random nature of the particle release feature.

Click

Global Evaluation 2

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Locate the section. In the table, enter the following settings:


Expression	Unit	Description
abs((Iref-cpt.pcnt1.It)/Iref)*100	1	Relative error (%) for the transmitted current
abs((yref-minop1(V,y))/yref)*100	1	Relative error (%) for the location of the potential minimum
abs((Vref-minop1(V))/Vref)*100	1	Relative error (%) for the potential minimum

The relative errors should be comparable in magnitude to the column in Table 1.

Click