of 17
View
Online
Export
Citation
CrossMark
RESEARCH ARTICLE
|
JULY 29 2021
Critical implications of ion-surface energy accommodation
and neutralization mechanism in hollow cathode physics
Special Collection:
Physics of Electric Propulsion
Pablo Guerrero
;
Ioannis G. Mikellides
;
James E. Polk
;
Rosa Carmina Monreal
;
Daniel I. Meiron
J. Appl. Phys.
130, 043306 (2021)
https://doi.org/10.1063/5.0055824
Articles Y
ou May Be Interested In
Dynamic thermal behavior of polycrystalline LaB
6
hollow cathodes
J. Appl. Phys.
(August 2021)
Evidence of nonclassical plasma transport in hollow cathodes for electric propulsion
J. Appl. Phys.
(March 2007)
The discharge plasma in ion engine neutralizers: Numerical simulations and comparisons with laboratory
data
J. Appl. Phys.
(December 2010)
04 October 2023 17:59:10
Critical implications of ion-surface energy
accommodation and neutralization mechanism
in hollow cathode physics
Cite as: J. Appl. Phys.
130
, 043306 (2021);
doi: 10.1063/5.0055824
View Online
Export Citation
CrossMar
k
Submitted: 3 May 2021 · Accepted: 6 July 2021 ·
Published Online: 29 July 2021
Pablo Guerrero,
1,2
,
a)
Ioannis G. Mikellides,
3
James E. Polk,
3
Rosa Carmina Monreal,
4
and Daniel I. Meiron
1
AFFILIATIONS
1
Department of Aerospace (GALCIT), California Institute of Technology, Pasadena, California 91001, USA
2
UCLA Samueli School Of Engineering, University of California, Los Angeles, California 90095, USA
3
Electric propulsion group, Jet Propulsion Laboratory, Pasadena, California 91109, USA
4
Departamento de Física Teórica de la Materia Condensada, Universidad Autónoma de Madrid and IFIMAC, Cantoblanco, Madrid
E28049, Spain
Note:
This paper is part of the Special Topic on Physics of Electric Propulsion.
a)
Author to whom correspondence should be addressed:
p.guerrero.eng@gmail.com
ABSTRACT
Self-heating thermionic hollow cathodes are essential components in modern plasma thrusters. To fully understand their operation, three inter-
dependent physical domains must be considered: plasma discharge physics, thermal response of the cathode structure, and chemical evolution
of plasma exposed surfaces. In this work, we develop the first self-consistently coupled plasma
thermal
chemical simulation platform for
hollow cathode operation using lanthanum hexaboride (LaB
6
) and Xe and study its performance against our experimentally determined tem-
perature measurements. Results show that the customary assumptions of single-step resonant neutralization and full energy accommodation in
ion-surface collisions fail to reproduce our empirical observations. We propose a two-step neutralization mechanism that consists of resonant
neutralization to the first excited state of xenon followed by Auger de-excitation to the ground state, along with system specific accommodation
factors. In this way, the agreement between the results of the simulations and experiments was achieved. These fundamental processes could
govern neutralization in other cathode technologies where low work function emitters are employed and should therefore be accounted for in
physical models. In addition, the new simulation platform allows us to better estimate the equilibrium work function of LaB
6
hollow cathode
emitters. In the cathode studied here, we found that the effective work function is 2.25 eV, which is significantly lower than previous estimates,
and leads to better than expected cathode material performance with important implications for space missions.
Published under an exclusive license by AIP Publishing.
https://doi.org/10.1063/5.0055824
I. INTRODUCTION
Reaching further, going faster, and delivering heavier payloads
are the ultimate design goals of any space vehicle. In order to
improve the current technical capabilities and enable more ambi-
tious deep space missions, NASA is investing in high power electric
propulsion systems. Scalable arrays of Hall and ion thrusters are at
the heart of these propulsion systems.
Hollow cathodes are among the basic components necessary
for Hall and ion thruster operation, and therefore cathode lifetime
is a major determinant of thruster lifetime. Hollow cathodes rely
on low work function materials to abundantly emit electrons into
the thruster plasma, where they are used for the critical tasks of
plasma generation and neutralization of positively charged exhaust
ion beams. Understanding hollow cathode physics is of great
importance given their critical role in thruster operation, however,
that involves knowledge in three domains of physics, i.e., plasma,
thermal, and chemical. Additionally, to fully capture the intricacies
of such a coupled environment, the key role of plasma
surface
interactions must be addressed.
Thermionic emission of electrons is the most fundamental
mechanism regulating cathode discharges. Thermionic emission
depends on the work function of the emitting surface of the
cathode (i.e., the emitter or
insert,
as shown in
Fig. 10
), which in
Journal of
Applied Physics
ARTICLE
scitation.org/journal/jap
J. Appl. Phys.
130,
043306 (2021); doi: 10.1063/5.0055824
130,
043306-1
Published under an exclusive license by AIP Publishing
04 October 2023 17:59:10
turn is defined by the chemical state and morphology of that
surface. Hollow cathode thermionic emission is modeled with the
Richardson
Dushman law,
1
j
ther
(
f
,
T
)
¼
DT
2
exp

e
(
f

f
Schottky
)
k
B
T

,
(1)
This law relates the emission current density of electrons
j
ther
to the temperature at the emitting surface
T
and the
(temperature-independent) work function of the material
f
.
D
is a
material-dependent constant that accounts for the work function
temperature dependence,
k
B
is the Boltzmann constant,
f
Schottky
is
the reduction in the work function due to the external electric field
applied to the surface (the Schottky effect), and
e
is the electron
charge.
Low work function is essential for electron emitters because it
enables operation at a lower temperature, which reduces the evapo-
ration of the emitting material and thus increases their useful life-
time. Currently, two mature technologies are widely used in plasma
thrusters: barium dispenser cathodes and lanthanum hexaboride
cathodes (LaB
6
).
2
These two different technologies are capable of
significant electron emission at relatively low temperatures and they
both depend on the cathode plasma acting as the heat source to
maintain the emitter temperature. Barium dispenser cathodes rely
on a supply of a barium-containing compound in the pores of a
porous tungsten emitter to maintain a layer of adsorbed barium
and oxygen atoms on the emitting surface which lowers the tung-
sten work function. Lanthanum hexaboride hollow cathodes
employ polycrystalline tubes machined from rods of press-sintered
LaB
6
powder. The bulk polycrystalline material serves as the elec-
tron emitter, so no surface complex of adsorbed atoms is required
to maintain a low work function and the lifetime is limited by
evaporation of the bulk material, not the supply of an activating
substance. LaB
6
cathodes are also more resistant to oxidizing con-
taminants than tungsten cathodes, which can greatly reduce propel-
lant feed system purification and handling requirements. We have
used LaB
6
emitters in this study but the results of this work are
also applicable to other cathode technologies.
Because evaporation of the insert material or depletion of the
barium supply is one of the main lifetime limiters for hollow cath-
odes, knowledge of the operating temperature is critical. Lifetime is
extremely sensitive to operating temperature because of the expo-
nential dependence of the evaporation rate on temperature. For
example, a 100

C higher insert temperature in a LaB
6
cathode can
result in a factor of ten reduction in lifetime.
3
Hollow cathodes are
self-heating devices that rely on their internal plasma discharge to
sustain the temperature required for emission. As a result, the
insert temperature is not directly controlled and the processes
involved in establishing the equilibrium temperature must be
understood to estimate the lifetime of these devices.
The operating temperature of hollow cathodes is established
by the interrelation between the internal plasma discharge distribu-
tion, the cathode structure thermal behavior, and the chemical state
of the emitter surface. Additionally, plasma
material interactions
determine heat fluxes and chemical evolution of the insert which
must be considered as well. The internal plasma discharge
characteristics can only be properly defined when the net electron
emission of the insert is known. This requires an accurate defini-
tion of the work function and temperature along the insert, in addi-
tion to accounting for the Schottky effect, space charge limitations,
and the return currents of electrons and ions to the cathode inter-
nal surfaces.
Hollow cathode modeling is typically simplified by neglecting
the thermal response of the cathode altogether and substituting
that part of the physics with measured temperatures and assumed
constant values for the insert work function.
4
7
Additionally, few
researchers
8
11
have coupled an extensive hollow cathode thermal
model with global plasma simulations, in 2D, so the true thermal
response of the device during self-heating operation and its depen-
dencies on the plasma loads have remained elusive.
In the present work, we address the limitations of previous
studies by first creating a high fidelity thermal model of the
cathode, second, incorporating a new approach to elucidate the
effective work function, and third, introducing limitations in the
energy accommodation and employing an alternative neutralization
modeling approach.
We introduce a novel high fidelity cathode thermal model and
validate it using a heat source governed by simple and well known
physics (Joule heating, see
Appendix A
). This provides two main
advantages: first, improved model accuracy, and second, the ability to
validate the thermal model independent of the plasma heating pro-
cesses which we are studying. Specifically, in our simulation plat-
form, we have employed the state-of-art, 2D axisymmetric, Orificed
Cathode (OrCa2D) plasma solver
7
,
12
14
to compute the plasma
structure, which provides a quantification of the net return current,
Schottky effect, and space charge limitations. We successfully
coupled the plasma and thermal solvers, providing a self-consistent
solution for the temperature distribution along the cathode insert.
This approach yields an accurate temperature profile along the insert,
which is then used by the plasma solver to compute the plasma
structure. The final temperature profile allows us to better estimate
the uneven evaporation of insert material and therefore to infer the
insert geometry evolution and the time to failure.
Furthermore, we describe a novel approach to include the
work function in our model. In our previous work,
15
we learned
that the LaB
6
hollow cathode work function at steady state cannot
be assumed to be that measured for vacuum diode cathodes,
2.67 eV.
1
,
16
A preliminary analysis suggested that the effective work
function for the cathode used here is

2
:
19 eV.
15
In this work, we
have used empirically determined plasma properties to precisely
quantify the effective work function of the material. In addition,
motivated by our previous findings, we studied the effect of work
function distributions along the insert on the self-consistent
solution.
Last and most importantly, we discovered that the commonly
assumed full energy accommodation of ion
surface interactions
and single-step resonant neutralization
17
yield results that disagree
with empirical temperature measurements. We introduce an alter-
native two-step neutralization mechanism and an estimation for
the energy accommodation threshold in our system. These alterna-
tive models reproduce very closely the empirical observations, and
we recommend that these phenomena be included in further work
on hollow cathode models.
Journal of
Applied Physics
ARTICLE
scitation.org/journal/jap
J. Appl. Phys.
130,
043306 (2021); doi: 10.1063/5.0055824
130,
043306-2
Published under an exclusive license by AIP Publishing
04 October 2023 17:59:10
II. METHODS
A. Overview of the coupled self-consistent simulation
approach
The overall strategy used to couple the different simulation
algorithms used in this work is shown in
Fig. 1
. MATLAB
®
was
used as the main framework from which OrCa2D and COMSOL
simulation packages are controlled, and where the coordination
between models is established. The general process of finding a
plasma solution, which is consistent with the thermal characteris-
tics of the hollow cathode requires several global iteration steps.
An iterative coupling of the plasma and thermal models is
possible due to the slower time scale of the thermal response of the
cathode compared to the much faster plasma behavior. A global
iteration step starts when the plasma model is provided with an
initial guess for the temperature distribution along the insert and
the plasma solution is then obtained for the particular operating
point. At that point, the heat flux model computes the heat fluxes
from the return currents. Those fluxes are then input into the
thermal model where the updated temperature of the insert is
obtained. Then, MATLAB
®
compares the temperature provided to
the plasma model at the beginning of the global iteration step with
the temperature profile that the thermal model has generated using
the plasma heat fluxes. If they are not the same, the merging func-
tion
g
1
(
T
(
z
)) shown in
Fig. 1
, Eq.
(2)
creates a new temperature
profile to input into OrCa2D, which is closer to the one that
COMSOL computed in the previous global iteration step.
T
0
i
þ
1
(
z
)
¼
T
i
(
z
)
if
j
T
i
(
z
)

T
0
i
(
z
)
j
,
dT
cap
8
z
,
T
0
i
(
z
)
þ
T
i
(
z
)

T
0
i
(
z
)

=
Fif
j
T
i
(
z
)

T
0
i
(
z
)
j
dT
cap
8
z
,

(2)
where
T
¼
T
insert
,
dT
cap
¼
parameter, and
F
¼
max(
j
T
i
(
z
)

T
0
i
(
z
)
j
)/
dT
cap
.
In order to input
T
0
insert
(
z
) into OrCa2D, a fourth order poly-
nomial fit is used. Then, a new global iteration can be executed.
The code finishes when the temperature input into OrCa2D and
the one computed by the COMSOL heat transfer module differ by
less than a small tolerance. The root mean square error (RMSE) of
both profiles is computed and convergence is defined when this
value is less than 10

C. At that point, we conclude that the plasma
solution is consistent with the thermal model of the cathode. With
this self-consistent solution, we are able to fully resolve
T
insert
(
z
),
along with a high fidelity estimate of the net return current and
Schottky effect.
7
The details of the chemical evolution of the
surface are captured by the input of a work function profile along
the emitter surface and the definition of the heat fluxes that the
plasma produces and input in the thermal model.
B. Hollow cathode thermal model
Thermal modeling was performed using COMSOL
Multiphysics
®
V5.3a. The cathode geometry was built using
SOLIDWORKS
®
3D CAD 2016 and imported into COMSOL. Two
Physics
modules were used to model heat transfer by conduction
and radiation between components,
Heat Transfer with
Surface-to-Surface Radiation
and
Heat Transfer in Thin Shells.
The
meshing approach was based on a parameterization using a combina-
tion of maximum sized elements in the edges of the domains and
extremely fine refinement
for the interior of the domain. The
thermal models were written in batch mode in Java, enabling
MATLAB
®
to dynamically create and execute them with any configu-
ration (geometry, material properties, boundary conditions, mesh con-
figuration and refinement and modeling details: parameters used for
thermal contact, convergence criter
ia,etc.).Anexampleofthesolution
of this strategy is shown in
Fig. 2
.The
baseline
thermal model was
validated independently of the plasma. This model uses the total
hemispherical emissivity of tantalum (
ε
Ta
) with an additional 0.2 to
define the radiation shielding. For more details about thermal model-
ing, see
Appendix A
.
C. Hollow cathode plasma model
The physics models, conservation equations, and numerical
methods in OrCa2D have been described in detail in previous
articles
7
,
12
14
and will only be described briefly here. The code
solves the conservation laws for three species in the partially
ionized gas: electrons, xenon ions, and xenon neutrals. It is
assumed that only singly charged ions are present and that quasi-
neutrality prevails except inside sheaths which are handled with
appropriate boundary conditions. The Navier
Stokes equations are
solved for the neutral gas only inside the cathode up to a
transi-
tion boundary
at which the method to obtain the solution changes
to a collisionless approach that assumes neutrals follow straight-line
trajectories.
18
The Euler equations for mass and momentum of
ions are solved in the entire computational domain. A separate
energy equation is solved for the ions, allowing for distinct temper-
atures for the two heavy species. Ionization, charge exchange, and
electron
ion collisions are accounted for in the equations and
modeled as either source or drag terms. The solution for the elec-
trons is obtained from a combination of Ohm
s law, energy, and
current conservation equations. An idealized model for the anoma-
lous enhancement in the resistivity, which is now known to occur
in the cathode plume
19
,
20
due to the excitement of ion acoustic
FIG. 1.
Overview of the strategy used to compute the
coupled solution that includes plasma, thermal, and heat
flux models.
Journal of
Applied Physics
ARTICLE
scitation.org/journal/jap
J. Appl. Phys.
130,
043306 (2021); doi: 10.1063/5.0055824
130,
043306-3
Published under an exclusive license by AIP Publishing
04 October 2023 17:59:10
turbulence (IAT), is included in Ohm
s law. The model is based on
the formulations of Sagdeev and Galeev (S&G)
21
and is described
in more detail in a paper by Mikellides
et al.
22
For more details
about the plasma model, see
Appendix B
.
D. Heat flux model
The inner region plasma consists of three fluids
neutrals,
electrons, and ions. Depending on the operating point of the
cathode, different fluxes of those species reach the internal cathode
surfaces exposed to the plasma. These fluxes heat the cathode, thus
maintaining operating temperatures. The heat fluxes are computed
using the following equations:
_
q
ther
¼
j
ther
f
,
(3)
_
q
e
¼
j
e

return
(2
T
e
þ
f
),
(4)
_
q
i
¼
j
i

return
α
0
EAC
f
(
β
)(
f
sheath
þ
T
e
=
2)
þ
(1

α
ISC
)
E
NC

, (5)
where
_
q
ther
is a cooling heat flux term due to emitted electrons
going over the work function barrier and being replaced by elec-
trons at the Fermi level.
_
q
e
and
_
q
i
are heating terms associated with
the return current of electrons and ions.
j
ther
is the thermionically
emitted electron current density defined in Eq.
(1)
.
j
e

return
is the
return current density of electrons.
j
i

return
is the return current
density of ions.
T
e
is the electron temperature at the sheath edge.
f
sheath
is the potential of the sheath.
α
ISC
is the ion survival proba-
bility.
α
0
EAC
is the energy accommodation coefficient for
normal-to-the-surface angle of incidence
β
.
f
(
β
) is the function
that models the dependence of the energy accommodation coeffi-
cient with ion-surface incidence angle.
E
NC
is the energy associated
with the neutralization event, with a maximum possible value equal
to the ionization potential of the first excited state of the gas minus
the work function of the surface.
In the cathode used in this study, the main thermal load con-
sists of ions returning to the cathode, accelerated in the sheath and
pre-sheath. The main cooling load comes from the thermionic
emission of electrons into the plasma. The emitter and the orifice
plate receive approximately the same net heat load. The next rele-
vant heat load is received by the keeper. The electron heating effect
is insignificant given the low electron return current collected.
The total amount of energy that ions returning to the cathode
bring with them consists of the energy picked up as the ions are
accelerated through the potential of the sheath
j
i

return
f
sheath
, their
thermal energy
j
i

return
T
e
=
2 (due to the plasma potential drop in
the pre-sheath), and the ionization potential (
IP
) of the specific gas
used, in our case, xenon.
23
The exact amount of heat transferred to
a surface when an ion collides with it depends on the following
factors: ion survival probability, energy accommodation factor, and
neutralization mechanism.
The ion survival probability
α
ISC
quantifies the probability
that an ion will survive (not be neutralized) when it interacts with a
surface. In our case, the low energy of the noble gas ions returning
to the surfaces of the cathode (
,
50 eV) indicates that we can
assume 0% probability of ion survival rate (
α
ISC
¼
0).
24
The return current of electrons is estimated following classical
sheath theory.
23
Returning electrons heat the material by the work
function as they descend through the potential barrier of the
surface. They also transfer their kinetic energy, computed as the
flux weighted average kinetic energy.
10
1. Energy accommodation coefficient
The energy accommodation coefficient
α
EAC
depends on the
inelastic behavior of the interaction between ions and the atomic
structure of the surface; therefore, it is a function of the energy and
incident angle of the ions with respect to the normal of the surface
(
β
).
α
EAC
only affects the potential and thermal energy of the ions,
where
α
0
EAC
is the energy accommodation factor for normal ion
incidence and
f
(
β
) is the function that accounts for deviations with
FIG. 2.
Hollow cathode temperature distribution computed with the thermal model described in this work (example).
Journal of
Applied Physics
ARTICLE
scitation.org/journal/jap
J. Appl. Phys.
130,
043306 (2021); doi: 10.1063/5.0055824
130,
043306-4
Published under an exclusive license by AIP Publishing
04 October 2023 17:59:10
respect to the normal. In this work, we only use the canonical form
of
f
(
β
)
¼
cos(
β
) as the ions become merely perpendicular to the
surface due to the strong action of the sheath and, therefore, the
results are basically insensitive to the functional form of
f
(
β
).
The energy accommodation coefficient strongly depends on
the mass ratio (
μ
) between gas ions and surface atoms. The EAC is
different for ion
surface interaction and neutral
surface interac-
tion. To our knowledge, there are no data available in the literature
for the system Xe
þ
LaB
6
in the energy range of interest for this
work. In order to bound the effect of the EAC in our case, we
extrapolated the results of the only data available for the interaction
between low energy ions and polycrystalline surfaces.
25
In our case,
Xe ions interact with LaB
6
crystals. Depending on the crystal, the
interaction will occur with lanthanum atoms [when the crystal is
purely terminated in La, like the crystal (001)] or a molecular form
of LaB
6
, depending on how many boron atoms participate in the
interaction. In the case Xe
þ
LaB
6
,
μ
Xe
þ
LaB
6
¼
0
:
64 and in the
case Xe
þ
La,
μ
Xe
þ
La
¼
0
:
95. According to Shuvalov
s experimen-
tal data,
25
the EAC
Xe
þ
LaB
6
for ion energies lower than 50 eV
must be lower than 90% and higher than 75% (
Fig. 3
).
2. Neutralization model
The energy associated with the neutralization event is accounted
for in
E
NC
. There are two possible families of neutralization mecha-
nisms based on how many steps it takes for the neutralization event
to be completed: single or two steps. The selection of mechanism
depends on the surface internal electron density distribution where
the ions collide and on the energy level of ionized and first excited
neutral states of the incoming ion.
1. One-step neutralization model
. Among the single-step neu-
tralization mechanisms, only two are applicable to the case where
low energy ions impact solid surfaces: resonant tunneling neutrali-
zation and Auger neutralization. The leading theory describes reso-
nant neutralization as the dominant mechanism by which low
energy ions neutralize to the ground state of the neutral atom upon
interactions with conductive surfaces.
26
Auger neutralization can
still exist in this case, with a much smaller probability. Only when
the density of states (DOS) in the solid is not aligned with the
energy levels in the ions, does Auger neutralization takes a domi-
nant role. This is the case when He
þ
interacts with high work func-
tion metal surfaces. In our case, if the neutralization event between
Xe
þ
and LaB
6
occurs following a single step, it will most probably
do it by resonant neutralization. In this case, incident xenon ions
are neutralized as they come closer to the surface by tunneling an
electron from the appropriate energy band of the solid into an
energy hole of the ion, thus creating a ground state xenon neutral.
In this case, the energy transferred to the surface is the ionization
potential of xenon minus the work function of the surface
E
NC
¼
IP

WF
, where
IP
¼
E
(Xe
þ
)

E
(Xe
0
).
2. Two-step neutralization model
. The alternative neutraliza-
tion model consists of resonant neutralization to the first excited
state of xenon (Xe
*
) followed by Auger de-excitation to the ground
state of xenon (Xe).
27
This mechanism is only possible if the work
function of the material is low enough, thus allowing electrons
from the appropriate energy level in the material to fill the respec-
tive energy hole in the ion, creating an excited neutral. In the case
of xenon, the difference between the ionized and the first excited
state is
E
(Xe
þ
)

E
(Xe
*
)
¼
3
:
81 eV.
28
This means that only the
emitter surface can be modeled with this two-step neutralization
model, as the other surfaces are pure metals and their work func-
tion is greater than 3.81 eV. The probability of creating excited neu-
trals is higher than creating ground state atoms because electrons
have to tunnel through a lower potential barrier.
In the first step, the amount of heat exchanged between the ion
and the surface is due to the transfer of one electron per ion neutral-
ized,
E
(Xe
þ
)

E
(Xe
*
)

WF
. The second step is the so-called Auger
de-excitation process, in which the excess of energy contained in the
excited xenon neutral relaxes to the ground state xenon neutral by
transferring that energy to the solid. The electrons in the solid will
then be excited with this energy. Depending on the density of states
in the solid, a different population of electrons will be excited and a
portion of them will be emitted by means of Auger de-excitation.
We model that yield with a parameter,
γ
Auger
. In summary, the total
amount of energy deposited in the solid due to neutralization per
ion strike is
E
NC
¼
IP

WF

γ
Auger
E
(Xe
*
)

E
(Xe
0
)
ðÞ
.
As a side effect, electrons are emitted from the solid at a rate
j
A
e
¼
γ
Auger
j
i

return
. These electrons get accelerated away from the
emitter by the emitter sheath. This component was accounted for
in OrCa2D when computing the net charge exchange at the
emitter surface. To our knowledge, there are no calculations of
γ
Auger
nor has it been experimentally quantified it for the combina-
tion Xe
þ
LaB
6
. Values below
γ
Auger
¼
3 have been measured for
several noble gases at different charge states interacting with
tungsten.
29
III. RESULTS
We employed the novel scheme described above to study the
self-consistent solution of a LaB
6
cathode operating at
J
D
¼
25 A
and
_
m
Xe
¼
13 sccm with the experimentally measured values of
V
K
¼
4
:
2 V and
V
D
¼
25 V for the keeper and discharge voltage.
FIG. 3.
EAC as a function of ion to surface mass ratio and ion energy.
25
Journal of
Applied Physics
ARTICLE
scitation.org/journal/jap
J. Appl. Phys.
130,
043306 (2021); doi: 10.1063/5.0055824
130,
043306-5
Published under an exclusive license by AIP Publishing
04 October 2023 17:59:10
In the simulations presented here, the plasma convergence criterion
was set to 0.002% relative error.
A. Sub-optimal performance of the baseline
thermal model
The self-consistent solution
T
insert
(
z
)
i
¼
T
0
insert
(
z
)
i

is shown
in
Fig. 4
for the
baseline
thermal model. The modeling approach
for the heat fluxes assumes
α
0
EAC
¼
1 and
f
(
β
)
¼
1. The neutraliza-
tion mechanism is the single-step resonant neutralization to Xe
0
and the work function profile along the insert is constant. In this
case, the work function converged to 2.55 eV and the net return
current, net electron return current, and net ion return current to
the cathode were

20%,

5%, and

26% of the discharge current,
respectively. The net return current, net electron return current,
and net ion returned current to the emitter were

16%,

2%, and

18% of the discharge current, respectively.
Notably, the converged insert temperature profile does not
match the thermocouple measurements, with the solution predicting
higher temperature (

150

C) than that measured with thermocou-
ples. The thermocouple measurement approach was thoroughly ana-
lyzed in our previous work,
15
and we concluded that the discrepancy
observed in the coupled converged solution cannot be explained by
an inaccuracy of the temperature measurements. Therefore, we con-
clude that this discrepancy must be the result of modeling inaccura-
cies in the plasma, heat flux, or thermal models.
B. Sensitivity to thermal model assumptions
To understand the sensitivity of the thermal model to the dif-
ferent parameters that define it, we studied the impact of parameter
variations from the
baseline
thermal model. The most sensitive
parameter in the cathode thermal behavior is the value for the total
hemispherical emissivity of tantalum (
ε
Ta
).
ε
Ta
is used to model
the radiation properties of the cathode tube, the heater, and the
radiation shield.
ε
Ta
in the literature is obtained for pure polished
material. We observed darkening of different Ta cathode compo-
nents, indicating impurities, that may increase the actual value of
the emissivity. The darkening was particularly prominent on the
radiation shielding; therefore, we used a value of
ε
Ta
þ
0
:
2 for this
component in order to obtain the agreement shown in
Fig. 12
.
This value was also used in the
baseline
simulation shown in
Fig. 4
.
We created a new thermal model to understand how the
base-
line
thermal model validation is impacted by the increase in the
emissivity that defines the radiation properties for some of its com-
ponents. This new thermal model underestimates the global
thermal losses of the cathode, and therefore, for any given power
input, the temperature distribution of the inner region (insert
region) is colder than the experimental values, see
Fig. 5
. The only
difference between this thermal model and the
baseline
is that we
have added 0.4 to the
ε
Ta
from literature in the material properties
utilized to model the radiation shield, as opposed to 0.2 used in the
baseline
. The new thermal model contains artificial thermal losses
(higher emissivity than used in the
baseline
thermal model) and
therefore it is not physically realistic (inconsistent with validation
data obtained with heater-only experiment,
Fig. 5
).
In
Fig. 5
, we have simplified the validation results by removing
the absolute values of the experiments and simulations that we
have shown in
Fig. 12
. As TC1, TC2, and TC3 give very similar
results due to their physical proximity, we only show results for
TC2. We have also removed the two simulation series 100%
P
and
95%
P
, and only show the mean temperature simulated by both.
Results for this thermal model show that TC2 is cooler than the
experimental results by a minimum of

30

C and a maximum of

110

C. In addition, the new thermal model estimates

70

C
cooler temperatures for all power levels than those for the
baseline
.
In order to estimate the sensitivity of the coupled system to
inaccuracies in the thermal model, we varied two parameters and
ran the global simulation until the converged solutions of the
system were found. These included
ε
Ta
(20% increase with respect
to the literature values) for every Ta component, and
ε
Ta
used in
the definition of the radiation shield (0.2 increase above the value
used in the
baseline
). Results are shown in
Fig. 4
.
When we increase the value of the emissivity for modeling the
radiation shield by 0.4 above the literature values, we are effectively
overestimating the losses in the cathode as shown in
Fig. 5
by the
cooler response of the insert. However, that does not affect the
response of the coupled solution significantly (

6

C), see
Fig. 4
.
This occurs because in the coupled system, the overestimation of
the thermal losses creates a temperature solution that results in
lower thermionic emission. In this case, the plasma solution cannot
converge given that the total discharge current of 25 A cannot be
satisfied. Therefore, the new plasma solution (self-consistent with
this perturbed thermal model) is different from the one for the
baseline
thermal model.
A more significant overestimation of the losses was performed
by using 120% of the literature values for the Ta emissivity for
every Ta component in the cathode. The response of the converged
solution using this modified thermal model was still not signifi-
cantly cooler (

24

C). Therefore, the discrepancy between the
converged global simulation results and the experimental values
cannot be explained by any reasonable modification of the cathode
FIG. 4.
Converged temperature results for the
baseline
thermal model.
Sensitivity analysis results for EAC, thermal model,
V
K
, and neutralization
model. Experimental values shown for comparison.
Journal of
Applied Physics
ARTICLE
scitation.org/journal/jap
J. Appl. Phys.
130,
043306 (2021); doi: 10.1063/5.0055824
130,
043306-6
Published under an exclusive license by AIP Publishing
04 October 2023 17:59:10
thermal characteristics that result in an increase in the thermal
model heat losses.
C. Sensitivity to plasma model parameters
The discrepancy between the predicted and observed tempera-
ture values could also be related to inaccuracies in the plasma
solver results due to sensitivities in the input parameters defining
the cathode operating point. The plasma solution depends on four
experimentally measured variables and the modeling techniques.
The four empirical variables
_
m
Xe
,
J
D
,
V
D
, and
V
K
were mea-
sured with accuracies that do not affect the final solution signifi-
cantly. However, the solution of the coupled system strongly
depends on the value of
V
K
.
V
K
is ultimately used to find the work
function at every global iteration step. We studied the sensitivity to
V
K
by reducing it to 3.5 V (0
:
7 V less than the actual measure-
ment). Results are shown in
Fig. 4
.
The self-consistent insert temperature solution was reduced by

30

C in this case. This result suggests that inaccuracies in
V
K
modeling are not a likely explanation of the observed difference
between experimental and modeling results.
D. Sensitivity to work function distribution
In this work, we also studied the dependence of the
self-consistent solution on work function distributions along the
FIG. 5.
COMSOL thermal model validation results. Comparison between thermal model
baseline
and modified version with 0.4 addition to the emissivity of the radiation
shield.
FIG. 6.
Work function profiles used for the sensitivity analysis.
FIG. 7.
Converged temperature results for the
baseline
thermal model.
Sensitivity analysis for work function distributions. Experimental values shown
for comparison.
Journal of
Applied Physics
ARTICLE
scitation.org/journal/jap
J. Appl. Phys.
130,
043306 (2021); doi: 10.1063/5.0055824
130,
043306-7
Published under an exclusive license by AIP Publishing
04 October 2023 17:59:10
emitter [
f
(
z
)]. This study was motivated by the observation of dis-
coloration along the emitter, suggesting chemistry evolution due to
plasma interaction which could result in such axial variations in
the work function. The impact on the magnitude of the work func-
tion or its spatial dependence is difficult to capture experimentally,
as we cannot measure the work function or the chemical state of
the surface while the cathode is running. However, to get an initial
estimate of the sensitivity of the self-consistent solution to work
function distributions, we decided to study the effect of linear and
parabolic work function profiles whose mean magnitude is adjusted
according to the self-consistent model needs. The amplitude of
those profiles is 0.2 eV, value chosen arbitrarily to get an initial esti-
mation of the sensitivity of the solution to a plausible real scenario,
see
Fig. 6
. We use
α
0
EAC
¼
0
:
75 for all the work function profiles
studied.
The self-consistent solutions show that the temperature distri-
bution is not greatly affected by the work function profiles utilized
in this study. The net effect on the insert temperature for profile 1
in
Fig. 6
is a shifted temperature distribution

30

C cooler than
the solution with constant work function profile. Similarly for
profile 3, it is

25

C cooler. For profiles 2 and 4, the shift is
toward higher temperatures, by

20

C(
Fig. 7
).
E. Sensitivity to energy accommodation
The results also showed that the solution is indeed insensitive
to
f
(
β
) and if
α
0
EAC
¼
0
:
75 at every surface of the cathode, the con-
verged solution does not agree with the experimental values as
shown in
Fig. 4
. The temperature distribution at convergence is in
this case (

74

C) cooler than the one obtained for the
baseline
thermal model.
F. Simulations incorporating two-step neutralization
model are consistent with empirical data
A much better agreement was found when the two-step
neutralization model was employed (see Sec.
II
). In this case,
FIG. 8.
Evolution of the temperature input to OrCa2D (black curves) and output from COMSOL thermal model based on heat fluxes from the last iteration of OrCa2D
(colored curves).
Journal of
Applied Physics
ARTICLE
scitation.org/journal/jap
J. Appl. Phys.
130,
043306 (2021); doi: 10.1063/5.0055824
130,
043306-8
Published under an exclusive license by AIP Publishing
04 October 2023 17:59:10