# Tunable strong magnetic anisotropy in two-dimensional van der Waals
antiferromagnets

Feiping Xiao  School of Physics and Electronics, Hunan University, Changsha
410082, China  Qingjun Tong<EMAIL_ADDRESS>School of Physics and
Electronics, Hunan University, Changsha 410082, China

###### Abstract

We show that anisotropic energy of a 2D antiferromagnet is greatly enhanced
via stacking on a magnetic substrate layer, arising from the sublattice-
dependent interlayer magnetic interaction that defines an effective
anisotropic energy. Interestingly, this effective energy couples strongly with
the interlayer stacking order and the magnetic order of the substrate layer,
providing unique mechanical and magnetic means to control the
antiferromagnetic order. These two types of control methods affect distinctly
the sublattice magnetization dynamics, with a change of the ratio of
sublattice precession amplitudes in the former and its chirality in the later.
In moiré superlattices formed by a relative twist or strain between the
layers, the coupling with stacking order introduces a landscape of effective
anisotropic energy across the moiré, which can be utilized to create
nonuniform antiferromagnetic textures featuring periodically localized low-
energy magnons.

KEYWORDS: 2D magnets, antiferromagnetic order, van der Waals heterostructure,
moiré pattern

##  I blueIntroduction

The demand of down scale of information carriers has spurred tremendous
interest in two-dimensional (2D) layered magnets, because of their potential
applications in constructing atomically thin spintronic devices Burch2018 ;
Gong2019 ; Mak2019 ; Huang2020 ; Rahman2021 . Compared with 2D ferromagnets
Gong2017 ; Huang2017 ; Jiang2018 ; Huang2018 ; Wang2018 ; Deng2018 ; Song2018
; Klein2018 ; Wang2018m ; Cenker2021 , antiferromagnetic (AFM) materials are
more robust for information storage, owing to their immunity against external
magnetic disturbance Wang2016 ; Du2016 ; Lee2016 ; Kim2019 ; Kang20202 ;
Wang20212 ; Ni20212 ; Hwangbo2021 . In addition, AFM magnets possess two
degenerate elementary chargeless bosonic modes (magnons), operating at
terahertz regime and allowing for ultrafast transfer of spin information
Kampfrath2011 ; Cheng2014 ; Rezende2019 ; Xing2019 ; Li2020 . A prerequisite
for 2D AFM application is the existence of a strong and tunable magnetic
anisotropy Jungwirth2016 ; Jungwirth2018 . First, a strong magnetic anisotropy
is needed to stabilize 2D magnetism against thermal fluctuations Mermin1966 ,
and achieve terahertz AFM resonance frequency Keffer1952 . Second, a tunable
magnetic anisotropy can be utilized to manipulate AFM order, whose control is
usually quite difficult due to the absence of net magnetization. This
motivates the ongoing searching for 2D AFM magnets with tunable strong
magnetic anisotropy.

Figure 1: Schematics of the response of magnetic orders with weak uniaxial
anisotropy to different sublattice magnetic field. (a) A uniform sublattice
magnetic field stabilizes the FM order, while destabilizes the AFM order into
spin flop phase. (b) A staggered sublattice magnetic field pointing in
opposite direction stabilizes the AFM order with two degenerate AFM resonance
modes. (c) A general staggered sublattice magnetic field stabilizes the AFM
order with two nondegenerate AFM resonance modes.

Besides the material-determined intrinsic anisotropy, magnetic order in a
ferromagnet can also be stabilized and oriented by an external magnetic field,
achieving the spin polarized state (c.f. Figure 1a). In contrast, because of
its staggered nature, the AFM magnetic order can only be controlled by a
staggered sublattice magnetic field that points along opposite direction (c.f.
Figure 1b). This requires atomic-scale magnetic control, which, however, is
impracticable for conventional magnetic techniques. The recent advent of 2D
van der Waals layered magnets sheds new light on this intractable problem
Mak2019 ; Huang2020 . A unique advantage of van der Waals materials is the
appearance of an additional degree of freedom regarding interlayer atomic
registry, i.e. the vertical arrangement of atomic sites of adjacent layers
Geim2013 ; Tong2017 ; Seyler2019 ; Tran2019 ; Jin2019 ; Alexeev2019 ; Song2019
; Li2019 ; Chen2019 ; Enaldiev2020 ; Stern2021 ; Yasuda2021 ; Woods2021 ;
Kennes2021 ; Andrei2021 ; Wilson2021 . When the AFM sublattice sits
differently on a magnetic substrate layer, they would experience a distinct
interlayer magnetic interaction. In the Landau-Lifshitz theory Landau1935 ,
this sublattice dependent magnetic interaction defines an effective staggered
sublattice magnetic field that operating at atomic scale (c.f. Figure 1c).
This staggered effective field arising from the unique stacking degree of
freedom in van der Waals materials has not been considered in conventional
magnetic control schemes, such as exchange bias where AFM materials are mainly
used as passive pinning substrates Berkowitz1999 . This motivates us to
explore the possibility to control AFM order in van der Waals layered 2D
magnets in pursuit of tunable strong magnetic anisotropy.

Here, via stacking an AFM monolayer on a ferromagnetic (FM) substrate, we show
that the interlayer magnetic interaction introduces a staggered sublattice
magnetic field, which can be divided into an effective anisotropic field that
stabilizes the AFM order and an effective Zeeman field that breaks the
degeneracy of magnon modes. Interestingly, the effective anisotropic field
couples strongly with interlayer stacking order and magnetic order of the
substrate layer, and even can change sign that flips the AFM order. These two
types of couplings affect distinctly AFM magnetization dynamics, where an
interlayer shift changes the ratio of sublattice precession amplitudes while a
magnetic field changes its chirality. In the moiré superlattice formed by a
relative twist or strain between the magnetic layers, the stacking dependent
effective anisotropic field creates nonuniform magnetic textures, including
AFM skyrmions. The low-energy magnons feel a trapping potential contributed
from both the AFM texture and the landscape of the effective field, and are
periodically trapped in the moiré. Our findings, supported by first-principles
calculations on MnPS3/CrCl3 heterobilayer, are general in layered 2D magnets
(AFM MnPS3 homobilayer is studied in the discussion section).

##  II blueEffective magnetic anisotropy and its tunability

We first give a general analysis on how to create a tunable magnetic
anisotropy from the sublattice dependent interlayer magnetic interaction. The
model Hamiltonian for a commensurate bilayer consisting of an AFM top layer
and an FM bottom layer is

$\displaystyle H^{bi}=\sum_{\langle
i,j\rangle,l}(J^{l}\bm{m}^{l}_{i}\bm{m}^{l}_{j}-K^{l}_{0}m^{2}_{z,i})-\sum_{i,j}J^{\perp}_{i,j}(\bm{r})\bm{m}^{t}_{i}\bm{m}^{b}_{j}$
(1)

where $\bm{m}^{l}_{i}=\bm{M}^{l}_{i}/M^{l}_{0}$ is the normalized magnetic
moment at $i$-site of the $l$-th layer with $M^{l}_{0}$ being its magnitude
and $J^{l}$ ($J^{\perp}_{i,j}$) is the intralayer (interlayer) exchange
coupling. The interlayer part depends on the interlayer distance and inplane
shift $\bm{r}$ between the magnetic atoms of the two layers Tong2017 . The
summation of lattice sites includes the sublattice degree of freedom and runs
over nearest neighbor sites $\langle i,j\rangle$ in the intralayer part.
$K_{0}^{l}$ is anisotropy energy, which is assumed to be small for the AFM top
layer so that its magnetic order can be tuned.

In the Landau-Lifshitz theory, the magnetic moments of the top AFM layer at A
and B sublattice feel a local effective field
$\bm{H}^{eff}_{A/B}(\bm{r})=-\partial H^{bi}/\partial\bm{m}_{i,A/B}$.
Neglecting the exchange field that does not affect the magnetic anisotropy and
assuming the FM order in the bottom layer is uniform, i.e.
$\bm{m}^{b}_{j}=\bm{n}$, this field reads

$\displaystyle\bm{H}^{eff}_{A/B}(\bm{r})$ $\displaystyle=$
$\displaystyle\pm\tau_{0}K^{t}_{0}\bm{z}+[\pm\tau
K_{eff}(\bm{r})+Z_{eff}(\bm{r})]\bm{n},$ (2)

where $\tau_{0}=sgn[\bm{l}\cdot\bm{z}]$ and $\tau=sgn[\bm{l}\cdot\bm{n}]$ with
the AFM Néel vector $\bm{l}=(\bm{m}_{A}-\bm{m}_{B})/2$, and we have defined
$K_{eff}(\bm{r})=\sum_{j}[J^{\perp}_{j}(\bm{r}_{A})-J^{\perp}_{j}(\bm{r}_{B})]/2$
and
$Z_{eff}(\bm{r})=\sum_{j}[J^{\perp}_{j}(\bm{r}_{A})+J^{\perp}_{j}(\bm{r}_{B})]/2$.
From Eq. (2), one immediately recognizes that the discrepancy of the two
sublattice effective field $K_{eff}(\bm{r})$ acts as an effective anisotropic
energy that can stabilize the AFM order. While their average $Z_{eff}(\bm{r})$
acts as a uniform effective Zeeman energy that can break the degeneracy of the
AFM resonance modes and even lead to spin flop phase when sufficiently large.
Importantly, contrary to the original anisotropic term $K_{0}^{t}$ that is
fixed by the material, this effective anisotropic energy $K_{eff}(\bm{r})$,
depending on not only the orientation of the FM order $\bm{n}$ but also on the
relative distance $\bm{r}$ between the magnetic moments of the top and bottom
layers, is highly tunable. With the former being controlled by an applied
magnetic field and the latter by an interlayer relative shift Jiang20182 ,
these couplings provide unique possibilities to manipulate AFM order.

Figure 2: (a) Four magnetic configurations used to calculate the interlayer
interaction induced effective fields. (b) A shifting vector $\bm{r}$ that
characterizes interlayer atomic registry of a commensurate MnPS3/CrCl3
bilayer, which is defined in a unit cell with unit vectors $\bm{a}_{1/2}$. (c)
First-principles calculated effective anisotropic field $K_{eff}$ (blue dots)
and Zeeman field $Z_{eff}$ (red squares) in an AFM monolayer MnPS3 induced by
interlayer interaction from an FM monolayer CrCl3 for various stacking
registries along the long-diagonal in (b). Three high-symmetry configurations
are shown atop. The curves are plotted using Eqs. (S1) and (S2) in Supporting
Information 1, whose stacking dependence in a full unit cell are shown in (d)
and (e). $\bm{n}=\bm{z}$ is used in the CrCl3 layer.

The underlying physics of the interlayer interaction induced effective
anisotropy becomes more evident when one considers the AFM sublattice
precession dynamics, via generalizing directly Kittel’s original derivation in
the case of a uniform magnetic field to a staggered sublattice magnetic field
Keffer1952 . For an FM substrate with $\bm{n}=\bm{z}$, the AFM resonance
frequency and the ratio of the precession amplitudes between two sublattice
are (c.f. Supporting Information 2)

$\displaystyle\omega_{\pm}(\bm{r})$ $\displaystyle=$ $\displaystyle\gamma
Z_{eff}(\bm{r})\pm\gamma\\{[K_{0}+\tau K_{eff}(\bm{r})]$
$\displaystyle[2J+K_{0}+\tau K_{eff}(\bm{r})]\\}^{1/2},$
$\displaystyle\chi_{\pm}(\bm{r})$ $\displaystyle=$
$\displaystyle-1-[K_{0}+\tau K_{eff}(\bm{r})]/J\mp\\{[K_{0}+\tau
K_{eff}(\bm{r})]$ (3) $\displaystyle[2J+K_{0}+\tau
K_{eff}(\bm{r})]\\}^{1/2}/J,$

where $\gamma$ is the gyromagnetic ratio and the superscript $t$ is omitted.
These two ($\pm$) modes have opposite chiralities that the magnetization
precesses circularly in opposite manners. It is evident that the effective
anisotropic energy $K_{eff}(\bm{r})$ plays the same role as $K_{0}$ and can
enhance the resonance frequency. Interestingly, with a typical small $K_{0}$,
a stable AFM order requires $\tau K_{eff}(\bm{r})>0$, i.e. $K_{eff}(\bm{r})>0$
stables Néel order with $\tau=1$ and vice versa. The effective Zeeman field
$Z_{eff}(\bm{r})$ lifts the degeneracy of the two modes, and the low-energy
mode vanishes at a critical field value, beyond which spin flop phase
develops.

The effective fields arising from the interlayer magnetic interaction can be
quantitatively estimated from the energy discrepancy of certain magnetic
configurations (c.f. Figure 2a). In particular, the effective anisotropic
energy related to the discrepancy of the sublattice field can be calculated
from the energy difference of configurations with top AFM layers that have
opposite magnetic orders (I and II configurations in Figure 2a), which gives
$E_{II}-E_{I}=4K_{eff}$. On the other hand, the effective Zeeman energy
related to the average of the sublattice field can be measured by the one with
top FM layers (III and IV configurations), which gives
$E_{IV}-E_{III}=4Z_{eff}$.

In the following, we illustrate our results with a concrete example of a
magnetic heterobilayer composed of an AFM MnPS3 and an FM CrCl3, both of which
have been experimentally reported in its monolayer form and have almost
identical honeycomb lattice constants. MnPS3 has a very weak anisotropic
energy $K_{0}\sim 0.002J$ Okuda1986  and has been verified to be lack of long-
range AFM order in its monolayer form Ni2021 . CrCl3 crystal has weak in-plane
anisotropy and its magnetic order can be tuned to out-of-plane by a small
applied magnetic field, which is much weaker than the effective Zeeman field
at the high-symmetry stackings McGuire2017 . For a commensurate honeycomb
bilayer, the interlayer stacking order is characterized by a vector $\bm{r}$,
which measures the inplane shift of a top Mn atom to a bottom Cr atom defined
in a monolayer unit cell (c.f. Figure 2b).

Figure 2c shows the first-principles calculated effective anisotropic field
(blue dots) and Zeeman field (red squares) along the long-diagonal of the unit
cell, where there exist three high-symmetry configurations shown atop. At AB
and BA stackings, where one of the AFM sublattice sits on top of the FM
magnetic atom while the other is on top of the hollow center, the discrepancy
of the interlayer magnetic interaction at the two sublattices is largest.
Indeed, a strongest anisotropic field $K_{eff}$ is found, reaching as high as
$40K_{0}$. Meanwhile, their opposite signs reflect a strong coupling between
the anisotropic field and stacking order. At AA stacking, where the A and B
sublattices of the top and bottom layer align, the effective field from the
two sublattice matches, resulting in a zero effective anisotropic energy.
Different from the anisotropic field, the effective Zeeman field is always
positive that aligns with the substrate FM order, and approaches its maximum
at AA stacking. Because the interlayer distance is much larger than the
inplane atomic variation in a unit cell, one can well fit the datas using
analytic expressions with low harmonic functions (c.f. curves in Figure 2c and
see also Eqs. S(1) and S(2) in Supporting Information 1), whose variation in
the whole unit cell are shown in Figures 2d and e Tong2017 ; Bistritzer2011 .

Figure 3: (a) Stacking dependence of $z$-component of the Néel order
$l_{z}=(m^{z}_{A}-m^{z}_{B})/2$ (top) and AFM resonance frequency $\omega$
(bottom). The Néel order is flipped around the AB and BA regions, and vanishes
around the AA region that enters into spin flop phase. The effective fields
are adopted from Figure 2 with $\bm{n}=\bm{z}$. (b) Schematics of the
precession of the AFM sublattice magnetization and its magnetic orientation
for three high-symmetry configurations, which are related to each other via an
interlayer shift between the layers exerted by a mechanical force (labeled by
“F”). (c) A magnetic field (labeled by “B”) tuning the substrate FM order
changes the AFM order in the top layer and its dynamics.

The coupling between the effective anisotropy field and interlayer stacking
order suggests a mechanical way to control the AFM order and its dynamics by
shifting relatively the two layers. The top panel of Figure 3a shows the
$z$-component of Néel order as a function of interlayer translation, when the
effective fields lie out of plane. Near the AB and BA stackings, a uniform AFM
magnetization with $|l_{z}|=1$ is developed but with an opposite sign, i.e.
the Néel order flips, as schematically shown in Figure 3b. Near AA stacking,
the strong Zeeman field together with weak anisotropy field destabilizes the
AFM phase and favors a spin flop one with $l_{z}=0$. The lower panel of Figure
3a shows the stacking dependence of AFM resonance frequency, which is largest
at AB and BA stackings, indicating their best magnetic stability. Because the
effective Zeeman field is always positive, the lower-energy mode is always
$\omega_{-}$ with a fixed chirality. At AB and BA stackings, although
degenerate in energy, these two modes have opposite sublattice precession
ratio because of the flipped Néel order (see Eq. (3)).

The coupling to the FM order of the substrate layer provides another magnetic
way to control the AFM order and its dynamics (c.f. Figure 3c). Different from
the mechanical control, which only flips the sign of the effective anisotropic
field, reversing the FM order would flip the sign of the effective Zeeman
field as well. In this case, the lower-energy mode changes from $\omega_{-}$
to $\omega_{+}$, which have opposite chiralities.

Figure 4: (a) Schematics of a moiré superlattice with unit vectors
$\bm{b}_{1,2}$ formed by twisting a commensurate magnetic bilayer. (b) Energy
per moiré cell as a function of moiré periodicity for AFM MB and SK textures,
which are schematically shown in (c) and (d). The dashed line is the energy
for uniform AFM state. The effective fields are adopted from Figure 2 with
$\bm{n}=\bm{z}$. (e) Magnon miniband of an MB texture with $b=42a$. (f)
Wavefunction distribution of the lowest-energy magnon (black dot in (e)),
which is localized around the AA region of a moiré.

##  III blueAFM textures in magnetic moiré superlattices

The stacking dependent effective anisotropic field provides a unique
possibility to engineer AFM textures Gomonay2018  in moiré superlattice formed
by a relative twist or strain between the commensurate heterobilayer (c.f.
Figure 4a). The periodicity $b$ of a moiré is inversely proportional to the
misaligned angle and lattice mismatch, which can be further controlled by a
twist or strain Ribeiro2018 ; Bai2020 . Compared with nonmagnetic one Geim2013
; Tong2017 ; Seyler2019 ; Tran2019 ; Jin2019 ; Alexeev2019 ; Enaldiev2020 ;
Stern2021 ; Yasuda2021 ; Woods2021 , a magnetic moiré, composed of lattices
assigned with magnetic atoms, can be regarded as a superlattice of magnetic
moments Tong2018 ; Sivadas20182 ; Hejazi2020 ; Akram2021 ; Song2021 ; Xu2021 ;
Xie2021 . The smoothly changing atomic registry in a moiré then introduces a
landscape of effective anisotropic and Zeeman fields for the top AFM layer. We
assume that the superlattice is formed by two rigid lattices so that there
exists a linear mapping between local region in the moiré and registry-
dependent commensurate stacking defined in a monolayer unit cell. We notice
that when the substrate FM layer has a nonuniform magnetic texture, the
effective field would also show a spatial dependence even for a commensurate
atomic bilayer. This FM texture is possible when a weak out-of-plane magnetic
field is applied to tune the initially in-plane magnetic order in the CrCl3
layer.

In terms of a spatially varying field, the AFM monolayer in the magnetic moiré
can be modeled by an AFM monolayer imposed with laterally modulated
anisotropic and Zeeman fields,

$\displaystyle H=\sum_{\langle
i,j\rangle}J\bm{m}_{i}\bm{m}_{j}-[K_{0}+K_{eff}(\bm{R}_{i})]m^{2}_{z,i}-Z_{eff}(\bm{R}_{i})m_{z,i},$

where $\bm{R}_{i}$ is the atomic position in a moiré and we have assumed that
the magnetization in the substrate layer is along $\bm{z}$-direction. The
stable magnetic textures in a moiré can be obtained by relaxing from different
initial magnetic configurations using Landau-Lifshitz-Gilbert equation
Landau1935 ; Gilbert2004 . For a long-wavelength moiré, we find that the
uniform AFM magnetization becomes unstable. Instead, magnetic textures, like
AFM magnetic bubble (MB) and skyrmion (SK), become stable low-energy
configurations (schematically shown in Figures 4c and d). With the increase of
moiré periodicity, these AFM textures become more stable, as evidenced by the
decrease in their energy (c.f. Figure 4b). Because SK has an inplane vortex
structure, it costs more intralayer exchange energy and thus is energetically
higher than MB. Interestingly, because these AFM textures are generated and
stabilized by a landscape of the effective field, they can be driven by a
slight interlayer shift and regulated by an external magnetic field, which
provides new means to tune AFM textures (c.f. Supporting Information 3). These
moiré generated AFM textures avoid the stringent requirement of
Dzyaloshinskii-Moriya interaction Zhang2017 , which can be utilized to
engineer topological phases of various quasiparticles Gobel2017 ; Diaz2019 ;
Diaz2021 .

Usually the spatial modulation of the magnetic order in magnetic textures
leads to low-energy magnons localized at domain walls Yu2021 . Instead, in
moiré magnetic textures, there exists an additional spatially modulated
effective field (c.f. Supporting Information 4). This field landscape tends to
localize low-energy magnons around the AA region, where the magnetization is
most unstable. Figure 4e shows the calculated magnon miniband, where a low-
energy band is separated from others. Its wavefunction is localized around AA
region of a moiré (c.f. Figure 4f). When thermal fluctuation is considered,
magnetic orders in these regions are more fragile than the ones in other
regions.

##  IV blueDiscussion and conclusions

The effective anisotropic field that stabilizes AFM order arising from the
interlayer magnetic interaction is a general result in van der Waals magnets.
A recent experimental study on van der Waals AFM MnPS3 has shown that long-
range magnetic order persists from bulk to bilayer, while absent in monolayer
Ni2021 . This is consistent with our result that interlayer magnetic
interaction creates an effective anisotropy that stabilizes long-range
magnetic order, while which is absent in the monolayer case without an
interlayer interaction. For bilayer MnPS3, our first-principles calculations
show that the interlayer interaction induced effective anisotropic energy can
reach as high as 12$K_{0}$ (c.f. Supporting Information 5).

In summary, we have revealed a tunable strong magnetic anisotropy in van der
Waals magnets that can stabilize and manipulate 2D AFM order. This magnetic
anisotropy arises from a sublattice dependent interlayer magnetic interaction
that strongly couples with the interlayer stacking order and magnetic order of
the substrate layer. These couplings provide unprecedented mechanical and
magnetic controls over AFM order as well as its dynamics. In the moiré of a
magnetic heterostructure, the stacking dependent anisotropic field can be
utilized to engineer nonuniform AFM textures, in which low-energy magnons are
periodically trapped in the moiré. Our results point to a new route to explore
and manipulate 2D AFM magnetism, with potential applications in constructing
atomic-scale AFM spintronic devices and studying magnon condensation physics.

##  V blueSupporting Information

The Supporting Information includes details of the first-principles
calculations, AFM magnetization dynamics, mechanical and magnetic control over
the moiré AFM textures, localized magnons in moiré AFM textures, and effective
anisotropic field in bilayer MnPS3.

##  VI blueAcknowledgments

We thank Wang Yao for stimulating discussions at initial stage of this work.
This work is supported by the National Natural Science Foundation of China
(11904095), the National Key Research and Development Program of Ministry of
Science and Technology (2021YFA1200503), and the Fundamental Research Funds
for the Central Universities from China.

## References

  * (1) Burch, K. S.; Mandrus, D.; Park, J.-G. Magnetism in two-dimensional van der Waals materials. _Nature_ 2018, 563, 47-52.
  * (2) Gong, C.; Zhang, X. Two-dimensional magnetic crystals and emergent heterostructure devices. _Science_ 2019, 363, 706.
  * (3) Mak, K. F.; Shan, J.; Ralph, D. C. Probing and controlling magnetic states in 2D layered magnetic materials. _Nat. Rev. Phys._ 2019, 1, 646-661.
  * (4) Huang, B.; McGuire, M. A.; May, A. F.; Xiao, D.; Jarillo-Herrero, P.; Xu, X. Emergent phenomena and proximity effects in two-dimensional magnets and heterostructures. _Nat. Mater._ 2020, 19, 1276-1289.
  * (5) Rahman, S.; Torres, J. F.; Khan, A. R.; Lu, Y. Recent Developments in van der Waals Antiferromagnetic 2D Materials: Synthesis, Characterization, and Device Implementation. _ACS Nano_ 2021, 15, 17175-17213.
  * (6) Gong, C.; Li, L.; Li, Z.; Ji, H.; Stern, A.; Xia, Y.; Cao, T.; Bao, W.; Wang, C.; Wang, Y.; Qiu, Z. Q.; Cava, R. J.; Louie, S. G.; Xia, J.; Zhang, X. Discovery of intrinsic ferromagnetism in two-dimensional van der Waals crystals. _Nature_ 2017, 546, 265-269.
  * (7) Huang, B.; Clark, G.; Navarro-Moratalla, E.; Klein, D. R.; Cheng, R.; Seyler, K. L.; Zhong, D.; Schmidgall, E.; McGuire, M. A.; Cobden, D. H.; Yao, W.; Xiao, D.; Jarillo-Herrero, P.; Xu, X. Layer-dependent ferromagnetism in a van der Waals crystal down to the monolayer limit. _Nature_ 2017, 546, 270-273.
  * (8) Jiang, S.; Shan, J.; Mak, K. F. Electric-field switching of two-dimensional van der Waals magnets. _Nat. Mater._ 2018, 17, 406-410.
  * (9) Huang, B.; Clark, G.; Klein, D. R.; MacNeill, D.; Navarro-Moratalla, E.; Seyler, K. L.; Wilson, N.; McGuire, M. A.; Cobden, D. H.; Xiao, D.; Yao, W.; Jarillo-Herrero P.; Xu X. Electrical control of 2D magnetism in bilayer CrI3. _Nat. Nanotechnol._ 2018, 13, 544-548.
  * (10) Wang, Z.; Zhang, T.; Ding, M.; Dong, B.; Li, Y.; Chen, M.; Li, X.; Huang, J.; Wang, H.; Zhao, X.; Li, Y.; Li, D.; Jia, C.; Sun, L.; Guo, H.; Ye, Y.; Sun, D.; Chen, Y.; Yang, T.; Zhang, J.; Ono, S.; Han, Z.; Zhang, Z. Electric-field control of magnetism in a few-layered van der Waals ferromagnetic semiconductor. _Nat. Nanotechnol._ 2018, 13, 554-559.
  * (11) Deng, Y.; Yu, Y.; Song, Y.; Zhang, J.; Wang, N. Z.; Sun, Z.; Yi, Y.; Wu, Y. Z.; Wu, S.; Zhu, J.; Wang, J.; Chen, X. H.; Zhang, Y. Gate-tunable room-temperature ferromagnetism in two-dimensional Fe3GeTe2. _Nature_ 2018, 563, 94-99.
  * (12) Song, T.; Cai, X.; Tu, M. W.-Y.; Zhang, X.; Huang, B.; Wilson, N. P.; Seyler, K. L.; Zhu, L.; Taniguchi, T.; Watanabe, K.; McGuire, M. A.; Cobden, D. H.; Xiao, D.; Yao, W.; Xu, X. Giant tunneling magnetoresistance in spin-filter van der Waals heterostructures. _Science_ 2018, 360, 1214-1218.
  * (13) Klein, D. R.; MacNeill, D.; Lado, J. L.; Soriano, D.; Navarro-Moratalla, E.; Watanabe, K.; Taniguchi, T.; Manni, S.; Canfield, P.; Fernández-Rossier, J.; Jarillo-Herrero, P. Probing magnetism in 2D van der Waals crystalline insulators via electron tunneling. _Science_ 2018, 360, 1218-1222.
  * (14) Wang, Z.; Gutiérrez-Lezama, I.; Ubrig, N.; Kroner, M.; Gibertini, M.; Taniguchi, T.; Watanabe, K.; Imamŏglu, A.; Giannini, E.; Morpurgo, A. F. Very large tunneling magnetoresistance in layered magnetic semiconductor CrI3. _Nat. Commun._ 2018, 9, 2516.
  * (15) Cenker, J.; Huang, B.; Suri, N.; Thijssen, P.; Miller, A.; Song, T.; Taniguchi, T.; Watanabe, K.; McGuire, M. A.; Xiao, D.; Xu, X. Direct observation of two-dimensional magnons in atomically thin CrI3. _Nat. Phys._ 2021, 17, 20-25.
  * (16) Wang, X.; Du, K.; Liu, Y. Y. F.; Hu, P.; Zhang, J.; Zhang, Q.; Owen, M. H. S.; Lu, X.; Gan, C. K.; Sengupta, P.; Kloc, C.; Xiong, Q. Raman spectroscopy of atomically thin two-dimensional magnetic iron phosphorus trisulfide (FePS3) crystals. _2D Mater._ 2016, 3, 031009.
  * (17) Du, K.-Z.; Wang, X.-Z.; Liu, Y.; Hu, P.; Utama, M. I. B.; Gan, C. K.; Xiong, Q.; Kloc, C. Weak Van der Waals Stacking, Wide-Range Band Gap, and Raman Study on Ultrathin Layers of Metal Phosphorus Trichalcogenides. _ACS Nano_ 2016, 10, 1738-1743.
  * (18) Lee, J.-U.; Lee, S.; Ryoo, J. H.; Kang, S.; Kim, T. Y.; Kim, P.; Park, C.-H.; Park, J.-G.; Cheong, H. Ising-Type Magnetic Ordering in Atomically Thin FePS3. _Nano Lett._ 2016, 16, 7433-7438.
  * (19) Kim, K.; Lim, S. Y.; Kim, J.; Lee, J.-U.; Lee, S.; Kim, P.; Park, K.; Son, S.; Park, C.-H.; Park, J.-G.; Cheong, H. Antiferromagnetic ordering in van der Waals 2D magnetic material MnPS3 probed by Raman spectroscopy. _2D Mater._ 2019, 6, 041001.
  * (20) Kang, S.; Kim, K.; Kim, B. H.; Kim, J.; Sim, K. I.; Lee, J.-U.; Lee, S.; Park, K.; Yun, S.; Kim, T.; Nag, A.; Walters, A.; Garcia-Fernandez, M.; Li, J.; Chapon, L.; Zhou, K.-J.; Son, Y.-W.; Kim, J. H.; Cheong, H.; Park, J.-G. Coherent many-body exciton in van der Waals antiferromagnet NiPS3. _Nature_ 2020, 583, 785-789.
  * (21) Wang, X.; Cao, J.; Lu, Z.; Cohen, A.; Kitadai, H.; Li, T.; Tan, Q.; Wilson, M.; Lui, C. H.; Smirnov, D.; Sharifzadeh, S.; Ling, X. Spin-induced linear polarization of photoluminescence in antiferromagnetic van der Waals crystals. _Nat. Mater._ 2021, 20, 964-970.
  * (22) Ni, Z.; Haglund, A. V.; Wang, H.; Xu, B.; Bernhard, C.; Mandrus, D. G.; Qian, X.; Mele, E. J.; Kane, C. L.; Wu, L. Imaging the Néel vector switching in the monolayer antiferromagnet MnPSe3 with strain-controlled Ising order. _Nat. Nanotechnol._ 2021, 16, 782-787.
  * (23) Hwangbo, K.; Zhang, Q.; Jiang, Q.; Wang, Y.; Fonseca, J.; Wang, C.; Diederich, G. M.; Gamelin, D. R.; Xiao, D.; Chu, J.-H.; Yao, W.; Xu, X. Highly anisotropic excitons and multiple phonon bound states in a van der Waals antiferromagnetic insulator. _Nat. Nanotechnol._ 2021, 16, 655-660.
  * (24) Kampfrath, T.; Sell, A.; Klatt, G.; Pashkin, A.; M$\ddot{\mathrm{a}}$hrlein, S.; Dekorsy, T.; Wolf, M.; Fiebig, M.; Leitenstorfer, A.; Huber, R. Coherent terahertz control of antiferromagnetic spin waves. _Nat. Photon._ 2011, 5, 31-34.
  * (25) Cheng, R.; Xiao, J.; Niu, Q.; Brataas, A. Spin Pumping and Spin-Transfer Torques in Antiferromagnets. _Phys. Rev. Lett._ 2014, 113, 057601.
  * (26) Rezende, S. M.; Azevedo, A.; Rodríguez-Suárez, R. L. Introduction to antiferromagnetic magnons. _J. Appl. Phys._ 2019, 126, 151101.
  * (27) Xing, W.; Qiu, L.; Wang, X.; Yao, Y.; Ma, Y.; Cai, R.; Jia, S.; Xie, X. C.; Han, W. Magnon Transport in Quasi-Two-Dimensional van der Waals Antiferromagnets. _Phys. Rev. X_ 2019, 9, 011026.
  * (28) Li, J.; Wilson, C. B.; Cheng, R.; Lohmann, M.; Kavand, M.; Yuan, W.; Aldosary, M.; Agladze, N.; Wei, P.; Sherwin, M. S.; Shi, J. Spin current from sub-terahertz-generated antiferromagnetic magnons. _Nature_ 2020, 578, 70-74.
  * (29) Jungwirth, T.; Marti, X.; Wadley, P.; Wunderlich, J. Antiferromagnetic spintronics. _Nat. Nanotechnol._ 2016, 11, 231-241.
  * (30) Jungwirth, T.; Sinova, J.; Manchon, A.; Marti, X.; Wun-derlich, J.; Felser, C. The multiple directions of antiferromagnetic spintronics. _Nat. Phys._ 2018, 14, 200-203.
  * (31) Mermin, N. D.; Wagner, H. Absence of Ferromagnetism or Antiferromagnetism in One- or Two-Dimensional Isotropic Heisenberg Models. _Phys. Rev. Lett._ 1966, 17, 1133 .
  * (32) Kittel, C. Theory of Antiferromagnetic Resonance. _Phys. Rev._ 1951, 82, 565.
  * (33) Zhang, C.; Chuu, C.-P.; Ren, X.; Li, M.-Y.; Li, L.-J.; Jin, C.; Chou, M.-Y.; Shih, C.-K. Interlayer couplings, Moiré patterns, and 2D electronic superlattices in MoS2/WSe2 hetero-bilayers. _Sci. Adv._ 2017, 3, e1601459.
  * (34) Tong, Q.; Yu, H.; Zhu, Q.; Wang, Y.; Xu, X.; Yao, W. Topological mosaics in moiré superlattices of van der Waals heterobilayers. _Nat. Phys._ 2017, 13, 356-362.
  * (35) Seyler, K. L.; Rivera, P.; Yu, H.; Wilson, N. P.; Ray, E. L.; Mandrus, D. G.; Yan, J.; Yao, W.; Xu, X. Signatures of moiré-trapped valley excitons in MoSe2/WSe2 heterobilayers. _Nature_ 2019, 567, 66-70.
  * (36) Tran, K.; Moody, G.; Wu, F.; Lu, X.; Choi, J.; Kim, K.; Rai, A.; Sanchez, D. A.; Quan, J.; Singh, A.; Embley, J.; Zepeda, A.; Campbell, M.; Autry, T.; Taniguchi, T.; Watanabe, K.; Lu, N.; Banerjee, S. K.; Silverman, K. L.; Kim, S.; Tutuc, E.; Yang, L.; MacDonald, A. H.; Li, X. Evidence for moiré excitons in van der Waals heterostructures. _Nature_ 2019, 567, 71-75.
  * (37) Jin, C.; Regan, E. C.; Yan, A.; Utama, M. I. B.; Wang, D.; Zhao, S.; Qin, Y.; Yang, S.; Zheng, Z.; Shi, S.; Watanabe, K.; Taniguchi, T.; Tongay, S.; Zettl, A.; Wang, F. Observation of moiré excitons in WSe2/WS2 heterostructure superlattices. _Nature_ 2019, 567, 76-80.
  * (38) Alexeev, E. M.; Ruiz-Tijerina, D. A.; Danovich, M.; Hamer, M. J.; Terry, D. J.; Nayak, P. K.; Ahn, S.; Pak, S.; Lee, J.; Sohn, J. I.; Molas, M. R.; Koperski, M.; Watanabe, K.; Taniguchi, T.; Novoselov, K. S.; Gorbachev, R. V.; Shin, H. S.; Fal’ko, V. I.; Tartakovskii, A. I. Resonantly hybridized excitons in moiré superlattices in van der Waals heterostructures. _Nature_ 2019, 567, 81-86.
  * (39) Song, T.; Fei, Z.; Yankowitz, M.; Lin, Z.; Jiang, Q.; H-wangbo, K.; Zhang, Q.; Sun, B.; Taniguchi, T.; Watanabe, K.; McGuire, M. A.; Graf, D.; Cao, T.; Chu, J.-H.; Cobden, D. H.; Dean, C. R.; Xiao, D.; Xu, X. Switching 2D magnetic states via pressure tuning of layer stacking. _Nat. Mater._ 2019, 18, 1298-1302.
  * (40) Li, T.; Jiang, S.; Sivadas, N.; Wang, Z.; Xu, Y.; Weber, D.; Goldberger, J. E.; Watanabe, K.; Taniguchi, T.; Fennie, C. J.; Mak, K. F.; Shan, J. Pressure-controlled interlayer magnetism in atomically thin CrI3. _Nat. Mater._ 2019, 18, 1303-1308.
  * (41) Chen, W.; Sun, Z.; Wang, Z.; Gu, L.; Xu, X.; Wu, S.; Gao, C. Direct observation of van der Waals stacking dependent interlayer magnetism. _Science_ 2019, 366, 983-987.
  * (42) Enaldiev, V. V.; Z$\acute{\mathrm{o}}$lyomi, V.; Yelgel, C.; Magorrian, S. J.; Fal’ko, V. I. Stacking Domains and Dislocation Networks in Marginally Twisted Bilayers of Transition Metal Dichalcogenides. _Phys. Rev. Lett._ 2020, 124, 206101.
  * (43) Stern, M. V.; Waschitz, Y.; Cao, W.; Nevo, I.; Watanabe, K.; Taniguchi, T.; Sela, E.; Urbakh, M.; Hod, O.; Shalom, M. B. Interfacial ferroelectricity by van der Waals sliding. _Science_ 2021, 372, 1462-1466.
  * (44) Yasuda, K.; Wang, X.; Watanabe, K.; Taniguchi, T.; Jarillo-Herrero, P. Stacking-engineered ferroelectricity in bilayer boron nitride. _Science_ 2021, 372, 1458-1462.
  * (45) Woods, C. R.; Ares, P.; Nevison-Andrews, H.; Holwill, M. J.; Fabregas, R.; Guinea, F.; Geim, A. K.; Novoselov, K. S.; Walet, N. R.; Fumagalli L. Charge-polarized interfacial superlattices in marginally twisted hexagonal boron nitride. _Nat. Commun._ 2021, 12, 347.
  * (46) Kennes, D. M.; Claassen, M.; Xian, L.; Georges, A.; Millis, A. J.; Hone, J.; Dean, C. R.; Basov, D. N.; Pasupathy, A. N.; Rubio, A. Moiré heterostructures as a condensed-matter quantum simulator. _Nat. Phys._ 2021, 17, 155-163.
  * (47) Andrei, E. Y.; Efetov, D. K.; Jarillo-Herrero, P.; MacDonald, A. H.; Mak, K. F.; Senthil, T.; Tutuc, E.; Yazdani, A.; Young, A. F. The marvels of moiré materials. _Nat. Rev. Mater._ 2021, 6, 201-206.
  * (48) Wilson, N. P.; Yao, W.; Shan, J.; Xu, X. Excitons and emergent quantum phenomena in stacked 2D semiconductors. _Nature_ 2021, 599, 383-392.
  * (49) Landau, L. D.; Lifshitz, E. M. On the Theory of the Dispersion of Magnetic Permeability in Ferromagnetic Bodies. _Phys. Z. Sowjetunion_ 1953, 8, 153-164.
  * (50) Berkowitz, A. E.; Takano, K. Exchange anisotropy-a review. _J. Magn. Magn. Mater._ 1999, 200, 552-570.
  * (51) Jiang, L.; Wang, S.; Shi, Z.; Jin, C.; Utama, M. I. B.; Zhao, S.; Shen, Y.-R.; Gao, H.-J.; Zhang, G.; Wang, F. Manipulation of domain-wall solitons in bi- and trilayer graphene. _Nat. Nanotechnol._ 2018, 13, 204-208.
  * (52) Okuda, K.; Kurosawa, K.; Saito, S.; Honda, M.; Yu, Z.; Date, M. Magnetic Properties of Layered Compound MnPS3. _J. Phys. Soc. Jpn._ 1986, 55, 4456-4463.
  * (53) Ni, Z.; Zhang, H.; Hopper, D. A.; Haglund, A. V.; Huang, N.; Jariwala, D.; Bassett, L. C.; Mandrus, D. G.; Mele, E. J.; Kane, C. L.; Wu, L. Direct Imaging of Antiferromagnetic Domains and Anomalous Layer-Dependent Mirror Symmetry Breaking in Atomically Thin MnPS3. _Phys. Rev. Lett._ 2021, 127, 187201.
  * (54) McGuire, M. A.; Clark, G.; KC, S.; Chance, W. M.; Jellison, G. E.; Jr.; Cooper, V. R.; Xu, X.; Sales, B. C. Magnetic behavior and spin-lattice coupling in cleavable van der Waals layered CrCl3 crystals. _Phys. Rev. Mater._ 2017, 1, 014001.
  * (55) Jung, J.; Raoux, A.; Qiao, Z.; MacDonald, A. H. _Ab initio_ theory of moiré superlattice bands in layered two-dimensional materials. _Phys. Rev. B_ 2014, 89, 205414.
  * (56) Gomonay, O.; Baltz, V.; Brataas, A.; Tserkovnyak, Y. Antiferromagnetic spin textures and dynamics. _Nat. Phys._ 2018, 14, 213-216.
  * (57) Ribeiro-Palau, R.; Zhang, C.; Watanabe, K.; Taniguchi, T.; Hone, J.; Dean, C. R. Twistable electronics with dynamically rotatable heterostructures. _Science_ 2018, 361, 690-693.
  * (58) Bai, Y.; Zhou, L.; Wang, J.; Wu, W.; McGilly, L.; Halbertal, D.; Lo, C. F. B.; Liu, F.; Ardelean, J.; Rivera, P.; Finney, N. R.; Yang, X.-C.; Basov, D. N.; Yao, W.; Xu, X.; Hone, J.; Pasupathy, A. N.; Zhu, X.-Y. Excitons in strain-induced one-dimensional moiré potentials at transition metal dichalcogenide heterojunctions. _Nat. Mater._ 2020, 19, 1068-1073.
  * (59) Tong, Q.; Liu, F.; Xiao, J.; Yao, W. Skyrmions in the Moiré of van der Waals 2D Magnets. _Nano Lett._ 2018, 18, 7194-7199.
  * (60) Sivadas, N.; Okamoto, S.; Xu, X.; Fennie, C. J.; Xiao, D. Stacking-Dependent Magnetism in Bilayer CrI3. _Nano Lett._ 2018, 18, 7658-7664.
  * (61) Hejazi, K.; Luo, Z.-X.; Balents, L. Noncollinear phases in moiré magnets. _Proc. Natl. Acad. Sci. U.S.A._ 2020, 117, 10721-10726.
  * (62) Akram, M.; LaBollita, H.; Dey, D.; Kapeghian, J.; Erten, O.; Botana, A. S. Moiré skyrmions and chiral magnetic phases in twisted CrX3 (X=I, Br, and Cl) bilayers. _Nano Lett._ 2021, 21, 6633-6639.
  * (63) Song, T.; Sun, Q.-C.; Anderson, E.; Wang, C.; Qian, J.; Taniguchi, T.; Watanabe, K.; McGuire, M. A.; St$\ddot{\mathrm{o}}$hr, R.; Xiao, D.; Cao, T.; Wrachtrup, J.; Xu, X. Direct visualization of magnetic domains and moiré magnetism in twisted 2D magnets. _Science_ 2021, 374, 1140-1144.
  * (64) Xu, Y.; Ray, A.; Shao, Y.-T.; Jiang, S.; Lee, K.; Weber, D.; Goldberger, J. E.; Watanabe, K.; Taniguchi, T.; Muller, D. A.; Mak, K. F.; Shan, J. Coexisting ferromagnetic-antiferromagnetic state in twisted bilayer CrI3. _Nat. Nanotechnol._ 2022 17, 143-147.
  * (65) Xie, H.; Luo, X.; Ye, G.; Ye, Z.; Ge, H.; Sung, S. H.; Rennich, E.; Yan, S.; Fu, Y.; Tian, S.; Lei, H.; Hovden, R.; Sun, K.; He, R.; Zhao, L. Twist engineering of the two-dimensional magnetism in double bilayer chromium triiodide homostructures. _Nature Physics_ 2022 18, 30-36.
  * (66) Gilbert, T. L. A phenomenological theory of damping in ferromagnetic materials. _IEEE Trans. Magn._ 2004, 40, 3443-3449.
  * (67) Zhang, X.; Zhou, Y.; Ezawa, M. Antiferromagnetic skyrmion: stability, creation and manipulation. _Sci. Rep._ 2016, 6, 24795.
  * (68) Göbel, B.; Mook, A.; Henk, J.; Mertig, I. Antiferromagnetic skyrmion crystals: Generation, topological Hall, and topological spin Hall effect. _Phys. Rev. B_ 2017, 96, 060406(R).
  * (69) Díaz, S. A.; Klinovaja, J.; Loss, D. Topological magnons and edge states in antiferromagnetic skyrmion crystals. _Phys. Rev. Lett._ 2019, 122, 187203.
  * (70) Díaz, S. A.; Klinovaja, J.; Loss, D.; Hoffman, S. Majorana bound states induced by antiferromagnetic skyrmion textures. _Phys. Rev. B_ 2021, 104, 214501.
  * (71) Yu, H.; Xiao, J.; Schultheiss, H. Magnetic texture based magnonics. _Phys. Rep._ 2021, 905, 1-59.


# Confirmation of Transit-Time Limited Field Emission in Advanced Carbon
Materials with Fast Pattern Recognition Algorithm

Taha Y. Posos<EMAIL_ADDRESS>Department of Electrical and Computer
Engineering, Michigan State University, MI 48824, USA  Oksana Chubenko
Department of Physics, Arizona State University, AZ 85287, USA  Sergey V.
Baryshev<EMAIL_ADDRESS>Department of Electrical and Computer Engineering,
Michigan State University, MI 48824, USA

###### Abstract

An accurate estimation of the experimental field-emission area remains a great
challenge in vacuum electronics. The lack of convenient means, which can be
used to measure this parameter, creates a critical knowledge gap, making it
impossible to compare theory to experiment. In this work, a fast pattern
recognition algorithm was developed to complement a field emission microscopy,
together creating a methodology to obtain and analyze electron emission
micrographs in order to quantitatively estimate the field emission area. The
algorithm is easy to use and made available to the community as a freeware,
and therefore is described in detail. Three examples of dc emission are given
to demonstrate the applicability of this algorithm to determine spatial
distribution of emitters, calculate emission area, and finally obtain
experimental current density as a function of the electric field for two
technologically important field emitter materials, namely an
ultrananocrystalline diamond film and a carbon nanotube fiber. Unambiguous
results, demonstrating the current density saturation and once again proving
that conventional Fowler-Nordheim theory, its Murphy-Good extension, and the
vacuum space-charge effect fail to describe such behaviour, are presented and
discussed. We also show that the transit-time limited charge resupply captures
the current density saturation behaviour observed in experiments and provides
good quantitative agreement with experimental data for all cases studied in
this work.

##  I Introduction

Many studies[1, 2, 3, 4, 5, 6] have convincingly demonstrated that the
electron emission from a large surface area field emission cathode placed in a
macroscopic electric field $E$ is not uniform. The emission area is only a
small portion of the total surface area of a cathode since most of the
emission is confined to a small number of emission spots randomly distributed
over the cathode surface. Therefore, a proper and thorough estimation of the
apparent emission area and a number of emission locations is essential to
quantify field emitters in terms of the current density $j$ and its variation
with the applied electric field $E$. The importance of developing such
methodologies is twofold. First is practical: an actual emission area needs to
be known to compare cathode materials produced by various or varied syntheses.
Second is fundamental: only a properly established $j$–$E$ (and not $I$–$V$)
relationship can be used to define the validity range of a classical Fowler-
Nordheim (FN) emission and clarify the role of other mechanisms that could
cause deviation from the FN emission, i.e. cause non-conventional behavior
that is being observed across a large body of experimental work[1, 2, 3, 4, 5,
6, 7, 8, 9, 10].

A simple and convenient way to measure distribution of electron emission sites
is using a luminescence (or phosphor) screen, also known as a scintillator.
These screens emit light when interact with electrons. When such a screen is
used as an anode in a field emission experiment, it magnifies and projects an
emission pattern formed on the cathode surface under the external field force.
When captured by a camera, it creates a micrograph. An experimental system can
be designed such that $I$–$V$ curves can be taken synchronously with
micrographs. Then, the micrographs can be used to extract the apparent
emission area and obtain the $j-E$ relation. This work is motivated by the
lack of a guided micrograph processing for emission area calculations. Here,
we present a fast image processing algorithm and demonstrate its application
for a thorough analysis of field emission data. Results are presented for
micrographs obtained from ultra-nano-crystalline dimond (UNCD) films and
carbon nano-tubes (CNTs) under the applied dc field using screens made of Ce-
doped yttrium aluminum garnet (YAG:Ce), which produces a bright green
luminescence line at 550 nm. The proposed image processing algorithm is also
applicable to any other phosphor screen that is able to produce spatially
separated features that are bright enough.

The algorithm is realized on the MATLAB platform and takes an advantage of its
strength in processing arrayed data. When compared to an earlier version
implemented in Mathematica,[1] where a server was required to process
extensive sets of micrographs, our present implementation performs 10 times
faster on a personal laptop. The image processing code is an open-source
software: the first and future releases as well as examples and the user
manual can be found in our GitHub page (see Ref. 11). To make the best use of
it, a thorough description of the mathematical background and implementation
details are given in Section II. In Section III, we emphasize the importance
of the established image processing framework to capture non-conventional
field emission behaviour of semiconducting nanodiamond and carbon nanotube
materials. Semiconductors and semimetals have long been known to violate the
classical field emission law.[12, 13, 14, 15, 16, 17, 1] Our approach makes it
possible to establish quantitative relation between the applied electric field
and experimental field emission current density and therefore allows to verify
theory against experiment directly. Results and findings of this work are
summarized in Section IV.

##  II The image processing algorithm

###  II.1 Feature Extraction

As shown in the examples below (Figs. 1A, 4A, and 6A), field emission patterns
can vary a lot, with an additional challenge being a bright halo background,
which complicates image analysis. However, after a closer look at these images
one notices that the patterns consist of bright spots separated spatially
enough to distinguish them visually. So, the strategy is to detect the
brightest pixel within each emission spot, which is further referred to as a
local maximum (LM). Then the number of electron emission sites is equal to the
LM count and the emission area can be estimated by assembling certain neighbor
pixels around LMs.

Figure 1: A) A typical micrograph obtained from an UNCD film under the applied
dc field.[1] A 450$\times$450 px2 image represents a projection of spatial
distribution of electron emission sites onto a YAG anode placed 106 $\mu$m
from a 4.4 mm diameter cathode. B) Detected LMs (shown with blue plus signs)
overlaid with emission spots shown on the micrograph. C) Emission pixels
(shown in blue), representing the projected emission area, overlaid with
emission spots shown on the micrograph.

Prior to numerical analysis, an RGB micrograph is converted into a gray-scale
format and represented by a 2D matrix of pixels. The intensity of each pixel
is given by an integer between 0 and 255. When a typical micrograph (Fig. 1A)
is represented as a 3D plot (Fig. 2), it can be seen that emission spots
appear as Gaussian peaks atop a certain background. LMs are the brightest
pixels of each Gaussian peak and they are also well separated in space.
Therefore, two features must be known for each pixel in order to define LMs:
pixel value and the distance to the nearest brighter pixel. Because pixels are
represented by integers, occurrence of more than one LM for each peak is
possible. In order to prevent this, a small random noise between 0 and 0.1,
excluding 0 and 0.1, is added to the image so that there are no identical
pixels in the data array. This procedure does not perturb the image because
the original image can be retrieved any time by rounding the matrix.

A pixel $a$ with spatial coordinates $(i_{a},j_{a})$, where $i_{a}$ is the row
number and $j_{a}$ is the column number in a 2D array of digitized image, will
be represented by $(p_{a},d_{a})$ in a feature space, where $p_{a}$ is the
intensity feature and $d_{a}$ is the distance feature. $p_{a}$ is just a pixel
value (an integer, typically between 0 and 255), which can be simply extracted
from a data array (before adding the noise).

Figure 2: A 3D plot of the micrograph shown in Fig. 1A.

A fast method of extracting the distance feature is searching for a brighter
neighbor in a certain neighborhood of each pixel, which is called the search
region[18, 19]. The search region must be large enough to enclose an entire
Gaussian-like peak, but small enough to enclose not more than one entire
Gaussian peak. All images presented in this work were 450$\times$450 px2.
Although the size of emission spots varied, two standard deviations of each
Gauss peak corresponded to 10 pixels or less. At the same time, the distance
between peaks was 20 pixels or more. Therefore, a circular search region was
chosen: centered around each pixel, the search region radius was set to 10
pixels. It should be noted that for different images or image sizes, this
value may have to be adjusted accordingly. The Euclidean distance between any
two pixels $a^{\prime}$ at $(i_{a^{\prime}},j_{a^{\prime}})$ and $a$ at
$(i_{a},j_{a})$ is defined as:

$d_{a^{\prime}a}=\sqrt{\left(i_{a^{\prime}}-i_{a}\right)^{2}+\left(j_{a^{\prime}}-j_{a}\right)^{2}}.$
(1)

If $s_{a}$ is the search region for a pixel $a$ carrying value $p_{a}$, and a
pixel $b$ carrying value $p_{b}$ contained inside the $s_{a}$ is the next
closest pixel to the pixel $a$ such that $p_{b}>p_{a}$, then the distance
$d_{a}$ assigned to the pixel $a$ is

$d_{a}=\sqrt{\left(i_{a}-i_{b}\right)^{2}+\left(j_{a}-j_{b}\right)^{2}}.$ (2)

If no other brighter pixel was found in a search region, the brightest pixel
is called the maximum in search region (MISR). There is a distance property
that can be defined for a MISR as the distance to another closest MISR that is
brighter than the former. Say, pixel $A$ with value of $p_{A}$ is the MISR in
its own search region $s_{A}$, and pixel $B$ with value of $p_{B}$ is the MISR
in its own search region $s_{B}$, and $p_{B}>p_{A}$ while $p_{B}$ is a MISR
closest to pixel $A$, then the distance $d_{A}$ assigned to pixel $A$ reads

$d_{A}=\sqrt{\left(i_{A}-i_{B}\right)^{2}+\left(j_{A}-j_{B}\right)^{2}}.$ (3)

The brightest pixel across the entire image is the global maximum (GM). There
is a distance value that is assigned to the GM. It is the distance to closest
MISR (regardless of its value). Consider, pixel $GM$ with value of $p_{GM}$ is
the global maximum and pixel $A$ with value of $p_{A}$ is a MISR, closest to
$GM$, the distance assigned to the $GM$ is

$d_{GM}=\sqrt{\left(i_{GM}-i_{A}\right)^{2}+\left(j_{GM}-j_{A}\right)^{2}}.$
(4)

After all these distances are extracted and recorded, the random noise added
earlier gets removed. The code stores both $(p_{a},d_{a})$ and $(i_{a},j_{a})$
pairs for every pixel. Scatter plots of intensity–distance arrays, which
correspond to emission micrographs from UNCD (Fig. 1A) and CNTs (Fig. 4A), are
presented in Figs. 3A and 5A, respectively. Such plots are called $decision$
$plots$, also known as $plots$ $of$ $features$ in machine learning literature.

Figure 3: A) Decision plot of extracted features for the micrograph obtained
from an UNCD cathode operated under the dc field. Green, black, and brown
dashed regions show locations of LMs, background pixels, and emission pixels
in the decision plot, respectively. There is no overlapping between an LM
cluster and a uniform background. B) A black curve shows the Gaussian decision
boundary. The pixels shown in red are classified as LMs. Figure 4: A) A
typical micrograph obtained from a CNT fiber under dc field. B) Detected LMs
(shown with blue plus signs) overlaid with emission spots shown on the
micrograph. C) Emission pixels (shown in blue), representing the projected
emission area, overlaid with emission spots shown on the micrograph. Figure 5:
A) Decision plot of extracted features for the micrograph obtained from a CNT
fiber under the dc field. Green, black, and brown dashed regions show
locations of LMs, background pixels, and emission pixels, respectively. The
red dashed region shows overlap between LMs and a background cluster. B) A
black curve is the applied Gaussian decision boundary. Red points are detected
LMs. Notice that not all points above the boundary are red; this is because
some false LMs were filtered out.

###  II.2 Decision Boundary

Both the density of emission spots and their intensity increase with applied
electric field. Moreover, the background (often non-uniform) arises due to a
large number of bright glowing spots. As a result, an analysis of decision
plots becomes more complicated. For example, two peaks can be located very
close to each other, so the distance feature of a fainter peak becomes very
close to or lower than the set distance boundary. Or, let’s say, there are
fairly faint but distinct emission spots atop a bright background. In this
case, there is no clear boundary between the intensity of emission spots and
the intensity of background pixels. Moreover, some background protrusions can
be mistaken as LMs. Therefore, an effective method for estimating an
appropriate decision boundary is required to identify LMs on the decision
plot.

A supervised machine learning using labeled data cannot be used for this class
of problems because labeling hundreds of images with hundreds of emission
spots is not feasible. In addition, micrographs vary a lot for different
materials, geometries, or evolve with a high dynamic range throughout a single
experiment. An unsupervised machine learning scheme might be used. However, as
seen from Figs. 3A and 5A, no distinct clusters form on the decision plots.
Therefore, a simple and robust approach using a tunable decision boundary of a
certain form is attempted here.

For pixels that come from a continuous low background (intense and high-
gradient background case will be discussed later), the algorithm finds a
brighter pixel less than a few pixels away. Thus, these background pixels are
characterized by a small distance $d_{a}$ and a small intensity $p_{a}$. As a
rule, the background pixels outnumber the emission-spot pixels. Therefore,
most of the background pixels lie near origin and appear as crowded clusters
encircled by black dashed loops in Figs. 3A and 5A. Some pixels that belong to
emission spots are also characterized by small distances (less than the radius
of the search region), but most of them are brighter than the background and
their pixel values are large, up to 255. Therefore, these pixels form
stretched clusters at the bottom of decision plots as shown by brown dashed
regions in Figs. 3A and 5A.

LMs are bright pixels separated from brighter pixels by large distances
(usually larger than the search radius). They are shown by dashed green
regions in decision plots. If the background is low and uniform (as in the
case shown in Fig. 1A), LMs can be easily separated from background pixels
(see the green dashed region in a corresponding decision plot in Fig. 3A).
However, when a stronger background appears with a distinct gradient across
the image plane as shown in Fig. 4A, there are relatively faint emission
pixels that are barely noticeable on a bright background. In Fig. 5A, these
pixels are shown inside the red dashed region and are intermixed with
background pixels. To separate LMs from the background, we use a Gaussian
decision boundary given by

$f(p_{a})=A\cdot e^{-\frac{(p_{a}-\mu)^{2}}{2\sigma^{2}}}+k,$ (5)

where $A$ is the amplitude, $\mu$ is the mean, $\sigma$ is the standard
deviation, and $k$ is the offset. Functions, which define decision boundaries
for micrographs in Figs. 1A and 4A, are shown in Figs. 3B and 5B,
respectively. The decision rule is as follows. Consider a pixel $a$ carrying
the pair of values $(p_{a},d_{a})$. If $d_{a}>f(p_{a})$, then the pixel $a$ is
a LM. The parameters $A$, $\mu$, $\sigma$, and $k$ can be adjusted for each
data set. If the background and LM pixels are mixed as shown by the red dashed
region in Fig. 5A, many background pixels are detected as LMs. One way to fix
it is to increase $\sigma$ to exclude background pixels from the LM list. In a
case when the density of emission spots is high (emission spots, and thus LMs,
are spatially very close to each other), decreasing $k$ helps to identify some
missing LMs. The distance feature $d_{a}$ was never greater than 300 for
450$\times$450 px2 images. Thus, the parameter $A$ was kept at the value of
300 for all data sets presented. Mean $\mu$ was set to zero for all data sets.
Any false LM can be further filtered out by applying the surface-fitting
method given in Sec. II.3. Therefore, a crude adjustment of the decision
boundary parameters is enough.

Decision boundaries shown with black curves in Figs. 3B and 5B were applied to
Figs. 3A and 5A, respectively. Then, false LMs were filtered out by the
surface-fitting method. A final list of detected LMs are labeled with red data
points on the decisions plots in Figs. 3B and 5B as well as with blue crosses
in Figs. 1B and 4B. These examples illustrate nearly perfect LM detection.

###  II.3 Emission Area

As it was already mentioned, each emission spot on a micrograph (e.g., Fig.
1A) appears on a 3D plot (Fig. 2) as a quasi-symmetric Gaussian with its
center at a LM. Once a LM is detected, a 2D Gaussian function can be used to
fit an intensity peak[20, 21]. The fitting function, which returns the
estimated pixel value $p_{ij}^{e}$ for a pixel at the position $(i,j)$, is
given by

$p_{ij}^{e}=A\cdot e^{-\frac{(i-i_{LM})^{2}+(j-j_{LM})^{2}}{2\sigma^{2}}}+C,$
(6)

where $A$ is the amplitude, $\sigma$ is the standard deviation, $C$ is the
offset to manage the background level, and $(i_{LM},j_{LM})$ is the position
of a given LM. Fitting parameters $A$, $\sigma$, and $C$ have to be determined
for each LM. A 15$\times$15 px2 region was chosen to fit each peak centered at
the position $(i_{LM},j_{LM})$, so $i_{LM}-7\leq i\leq i_{LM}+7$ and
$j_{LM}-7\leq j\leq j_{LM}+7$. Such fitting region is optimal for our images
because it ensures that only one emission spot is included. The fitting
parameters for each LM can be calculated via minimizing the following sum by
the least squares regression method

$\sum_{i}\sum_{j}(p_{ij}-p_{ij}^{e})^{2},$ (7)

where $p_{ij}$ is the original pixel value, $i$ and $j$ are the indices
running over the fit region of each LM.

Obviously, most of the bright pixels will fall within one standard deviation
of the mean. To classify a pixel $a$, which is located at $(i_{a},j_{a})$
within a fit region centered at a LM at $(i_{LM},j_{LM})$, as the $emission$
$pixel$ that contributes toward the emission area calculation, the following
condition must be satisfied

$(i_{a}-i_{LM})^{2}+(j_{a}-j_{LM})^{2}<\sigma^{2}.$ (8)

Alone, this condition is not enough to classify a pixel as the emission pixel.
For example, consider the case when a background pixel is incorrectly detected
and listed as a LM (we call it a false LM in Sec. II.2). Gaussian fit of the
region around such a false LM will yield a large standard deviation $\sigma$,
but small amplitude $A$. This allows for filtering out fake LMs by applying
threshold values for $A$ and $\sigma$. LMs, which satisfy the condition

$(\sigma>\sigma_{th})\lor(A<A_{th}),$ (9)

are not real LMs and they are discarded from the master LM list. Thus, they do
not contribute to the local maxima count and the emission area calculation.

For the micrographs presented in this work, $\sigma_{th}=7$ and $A_{th}=10$
are chosen. For cases when spots are faint, $A_{th}$ should be decreased to
make the spots detected. If background noise is high, $\sigma_{th}$ should be
decreased to produce reliable results. In Figs. 1C and 4C, the blue pixels
identified as the emission pixels are overlaid with the emission spots shown
in Figs. 1A and 4A, respectively. If the physical area of the image is known
(we typically work with square images with a side length equal to the diameter
of the cathode), then the apparent emission area is given by

$\begin{split}\text{Emission area (mm${}^{2}$)}&=\text{Image area
(mm${}^{2}$)}\\\ &\quad\times\frac{\text{Number of emission
pixels}}{\text{Total number of pixels}}.\end{split}$ (10)

Both Figs. 1A and 4A have 450$\times$450=202,500 pixels. There are 1,929
emission pixels in Fig. 1C and 349 emission pixels in Fig. 4C. The displayed
portion of a YAG:Ce screen in Figs. 1A and 4A has the area of 4.4$\times$4.4
mm2 and 12.4$\times$12.4 mm2, respectively. Then applying Eq. 10, the apparent
emission area is 0.814 mm2 for UNCD (Fig. 1) and 0.264 mm2 for CNTs (Fig. 4),
while the total cathode surface area available for emission is 15.2 mm2 for
UNCD and 4.71 mm2 for CNTs.

###  II.4 Micrographs with Intense Background

Figure 6: A) A typical micrograph obtained from a CNT fiber under dc field
with glowing high-gradient background. The emission domains appear as bright
peaks on the glow. Exact source of the glow is unknown. B) Detected LMs (shown
with blue plus signs) overlaid with emission spots shown on the micrograph. C)
Emission pixels (shown in blue), representing the projected emission area,
overlaid with emission spots shown on the micrograph. Figure 7: A) Decision
plot of extracted features for the micrograph obtained from a CNT fiber under
the dc field with glowing high-gradient background. Unlike the previous cases,
the background does not form a cluster. Instead, it is distributed over a wide
pixel range shown in a black dashed region. Green dashed region consists of
well separated LMs. Red region includes both background pixels and LMs, so the
boundary should be drawn so that that region includes candidate LMs. Any false
LMs are to be filtered out by Gaussian surface fitting. B) A black curve shows
the applied decision boundary. All points above the curve are candidate LMs.
Red points label finalized LM list after false LMs are filtered out by the
surface fitting.

Micrographs with intense and non-uniform backgrounds are extremely challenging
to process and analyze. One such example is shown in Fig. 6A. In this
micrograph, the global background changes from dark at the top to very bright
at the bottom, where the detection of emission spots becomes less effective. A
corresponding decision plot is shown in Fig. 7A. Unlike previous cases (Figs.
3 and 5), the background pixels do not cluster near the origin, but instead
they are distributed over a wide range of pixel values (shown by a black
dashed region in Fig. 7A). There are well separated bright pixels,
corresponding to LMs, shown by a green dashed region. Because some of the
glowing background pixels have values comparable or higher than those of some
of the emission pixels, there is also a large region, shown by a red dashed
line, where LMs and background pixels are intermixed.

It is important to note that any LM, excluded during the decision-making
procedure, cannot be recovered at later stages of the image analysis. However,
any background pixel, identified and included as a LM, can be filtered out by
the surface-fitting procedure. Therefore, the decision boundary should be
$soft$, i.e. it should include all possible LM candidates (as opposed to a
$rigid$ boundary that would exclude as many background pixels as possible).
According to this approach, the decision boundary should be drawn on the
decision plot (Fig. 7) in such a way that the pixels from a red region are
considered as LM candidates. Any false LM will be likely filtered out. One
major drawback of applying the soft boundary is the increased computation time
due to increased number of Gaussian fitting steps required for increased
number of LMs. With the decision boundary applied in Fig. 7B, all pixels above
the black curve are LM candidates. Fig. 7B highlights in red a finalized array
of LMs (after most or all false LMs were removed). Same points are overlaid on
the original image as blue crosses in Fig. 6B. Nearly perfect agreement can be
seen upon visual inspection. Finally, Fig. 6C shows in blue the calculated
apparent emission area. One can see how effectively the entire pattern
recognition workflow performs even for images highly distorted by the
background. The calculated number of emission pixels is 533 out of 202,500
pixels. The field of view in Fig. 6A is 11.2$\times$11.2 mm2 yielding an
emission area of 0.33 mm2 while the total cathode area is 15.2 mm2.

##  III Results and Discussions

Figure 8: A) SEM image of the UNCD film grown on a Ni/Mo/SS substrate. B) SEM
image of the CNT fiber. The twisted and folded yarn structure is clear from
the image.

In this section, we use the described algorithm to extract the electron
emission area projected on the YAG screen, or the apparent emission area, and
estimate the experimental current density as a function of the electric field.
We analyze three sets of field emission micrographs. Each set contains tens of
micrographs taken with regular intervals as the applied voltage (and the
electric field) was swept up and down. The current was measured concurrently.
The first set is obtained from the $n$-type UNCD film grown on a planar
stainless steel (SS) substrate with a Ni/Mo buffer layer[1]. The second and
third sets are from the CNT folded rope taken before and after conditioning,
respectively. This sample was formed by braiding and twisting 90 $\mu$m
fibers, then folding the resulting yarn multiple times and inserting it in a
metal tube.[2] Scanning electron microscopy (SEM) images of the samples are
given in Fig. 8.

Electric field dependencies of the apparent emission area
$S_{\text{apparent}}$ computed for the considered data sets are presented with
blue dots in Fig. 9. In the same figure, red dots trace corresponding field
dependencies of the measured current $I_{\text{measured}}$. As it was already
mentioned above, the apparent emission area is not the actual emission area.
The apparent emission area estimated by the algorithm is the projected, or
magnified, emission area. The magnification occurs as a result of the non-zero
finite transverse momentum of emitted electrons as well as the transverse
component of the fringe electric field on the surface of the high-aspect-ratio
features. Therefore, the relation between the apparent emission area and the
real emission area can be obtained by introducing the magnification factor.
Although finding the factor is not a straightforward process (because it
depends on both internal material parameters and complicated surface
geometry), the lower and upper limits of the experimental current density can
be easily estimated.

The lower limit of the current density can be found by dividing the measured
current $I_{\text{measured}}$ by the total apparent emission area
$S_{\text{apparent}}$ estimated by the algorithm for a given electric field
$E$

$j_{\text{exp.min}}=\frac{I_{\text{measured}}}{S_{\text{apparent}}}.$ (11)

The lower current-density limits $j_{\text{exp.min}}$ for all three data sets
are shown with green dots in Fig. 10. It can be seen that the current density
$j_{\text{exp.min}}$ does not change significantly with the applied electric
field. This happens because the variation of the apparent emission area with
the applied electric field is almost identical to that of the measured current
(Fig. 9). Such behaviour is observed for all data sets and contradicts to the
Fowler-Nordheim law,[22] which predicts the exponential increase of the
current density with the applied electric field.

Further on, the upper limit of the current density can be calculated using an
estimated angle of the beam ray trajectory. If the angle is $\alpha$, then a
single emission site/domain with area $S_{\text{single}}$ appears on the
YAG:Ce screen as $S_{\begin{subarray}{c}\text{single}\\\
\text{apparent}\end{subarray}}$ given by

$S_{\begin{subarray}{c}\text{single}\\\
\text{apparent}\end{subarray}}=\pi(d\tan{\alpha})^{2},$ (12)

where $d$ is the inter-electrode gap. Geometrically, the angle $\alpha$ can be
obtained using the mean transverse energy (MTE) as[2]

$\tan{\alpha}=\frac{p_{x}}{p_{z}}=\sqrt{\frac{2\cdot\text{MTE}}{m_{0}c^{2}}}\frac{1}{v/c}\frac{1}{\gamma},$
(13)

where $m_{0}c^{2}=0.511\ \text{MeV}$ is the electron rest energy, $v/c$ is the
ratio of electron speed to speed of light, and $\gamma$ is the Lorentz factor.
At $1$ keV energy, $v/c$ and $\gamma$ are taken as $0.063$ and $1$,
respectively. MTE values are 100-200 meV[23, 24] for the UNCD; and for CNT its
Fermi energy of 50-100 meV[25] can be used. Then, the estimated values of the
angle is $0.9^{\circ}$ for the UNCD sample[26] and $0.4^{\circ}$ for the fiber
sample[2]. Then the inverse magnification factor is given by the ratio
$S_{\text{single}}/S_{\begin{subarray}{c}\text{single}\\\
\text{apparent}\end{subarray}}$ and the current density can be calculated as

$\begin{split}j_{\text{exp.max}}=\frac{I_{\text{measured}}}{\frac{S_{\text{single}}}{S_{\begin{subarray}{c}\text{single}\\\
\text{apparent}\end{subarray}}}\cdot S_{\text{apparent}}}.\end{split}$ (14)

Figure 9: The apparent emission area estimated using the algorithm (blue dots)
and the measured current (red dots) as a function of the applied electric
field A) for the UNCD film, B) for the CNT fiber before conditioning, and C)
for the CNT fiber after conditioning. Figure 10: The experimental current
densities (the lower limit is shown with green dots; the upper limit is shown
with blue dots) compared with the theoretical estimations (the elementary FN
equation is shown with magenta dots; the Murphy-Good extension is shown with
cyan dots) as well as the space-charge limited emission (dashed red line) and
the saturation-velocity limited emission (dashed orange line) as a function of
the applied electric field A) for the UNCD film, B) for the CNT fiber before
conditioning, and C) for the CNT fiber after conditioning.

To obtain the upper current density limit, the minimum possible size for
$S_{\text{single}}$ should be chosen in Eq. 14. Therefore, we assume
$S_{\text{single}}$ to be equal to the size of grain boundaries ($1$ nm by $1$
nm square[26]) and the size of a single CNT ($10$ nm in diameter circular
region) for the UNCD and the fiber, respectively. The resulting upper limits
for the three datasets are shown with blue dots in Fig. 10. A few important
conclusions can be made. First, experimental current densities are practically
constant in the studied range of electric fields. Therefore, an increase in
current with electric field can be explained by an increase in emission area
(see Fig. 9). Second, according to Dyke $et$ $al.$[27], above
$10^{7}\text{A}/\text{cm}^{2}$ current density level, the emission is limited
by the space charge, which causes decrease in effective field over the
cathode. The current density limit due to the space charge is shown by a red
dashed line for reference in Fig. 10. However, we find that the upper limit of
experimental current density is much lower than this value. Therefore, the
saturation behaviour cannot be explained by the space-charge effect.

We also compare our results with predictions given by the field emission
theory in both Fowler-Nordheim (FN) and Murphy-Good (MG) formulations. In FN
formulation, the image-charge effect is ignored and the potential barrier at
the material-vacuum interface can be approximated by a simple triangular
potential. In this case, the FN current density
$j_{\text{FN}}~{}(\frac{\text{A}}{\text{m}^{2}})$ is given by[22, 28]

$j_{\text{FN}}=A\frac{1}{\phi}(\beta E)^{2}e^{-\frac{B\phi^{3/2}}{\beta E}},$
(15)

where $A=1.54\times 10^{-6}\frac{\text{AeV}}{\text{V}^{2}}$, $B=6.83\times
10^{9}\frac{\text{V}}{\text{eV}^{3/2}\text{m}}$, $\beta$ is the unitless field
enhancement factor, $\phi~{}(\text{eV})$ is the material work function which
can be taken as 4.8 eV for both UNCD and fiber samples, and
$E~{}(\frac{\text{V}}{\text{m}})$ is the electric field. $\beta$ can be
obtained from experimental data by plotting it in FN coordinates where
$x$-axis is $1/E$ and $y$-axis is $\ln(j/E^{2})$. The plot forms a distinct
knee point for non-metal and semi-metal emitters[1, 2] that highlights the
saturation behaviour in a high-field region. After fitting a low-field region
with a line as in Ref. 2, the slope $m$ is given by

$m=-\frac{B\phi^{3/2}}{\beta},$ (16)

where $\beta$ can be extracted from. The resulting current density curves are
shown by magenta dots in Fig. 10.

A better theoretical approach is taking the image-force/charge effect into
account. In this case, the surface potential barrier is given by a rounded
triangular barrier and the current density is given by the MG current density
$j_{\text{MG}}~{}(\frac{\text{A}}{\text{m}^{2}})$[29, 30, 28, 31, 32, 33]

$j_{\text{MG}}=A\frac{1}{\text{t}_{\text{F}}^{2}\cdot\phi}(\beta
E)^{2}e^{-\text{v}_{\text{F}}\frac{B\phi^{3/2}}{\beta E}},$ (17)

where $A$ and $B$ are the same constants as before, and $\text{t}_{\text{F}}$
and $\text{v}_{\text{F}}$ are special field emission elliptic functions of
unitless Nordheim parameter $y$ given by[32, 30]

$y=k\cdot\frac{1}{\phi}\sqrt{\beta E},$ (18)

where $k=3.79\times
10^{-5}~{}\frac{\text{eV}\cdot\text{m}^{1/2}}{\text{V}^{\text{1/2}}}$, $E$ and
$\phi$ are in terms of $\text{V}/\text{m}$ and eV, respectively. Once it is
plotted in FN coordinates, it produces the slope $m$ given by

$m=-\text{s}_{\text{F}}\frac{B\phi^{3/2}}{\beta},$ (19)

where $\text{s}_{\text{F}}$ is the slope correction factor. Although
$\text{v}_{\text{F}}$, $\text{t}_{\text{F}}$, and $\text{s}_{\text{F}}$ are
not trivial to calculate, they had been tabulated well[34] and simple good
approximations are given by[35]

$\text{v}_{\text{F}}=1-y^{2}+\frac{1}{3}y^{2}\ln{y},$ (20)
$\text{t}_{\text{F}}=1+\frac{1}{9}(y^{2}-y^{2}\ln{y}),$ (21)
$\text{s}_{\text{F}}=1-\frac{1}{6}y^{2}.$ (22)

For a typical range of applied fields and the work function of 4.8 eV,
$\text{s}_{\text{F}}$ and $\text{t}_{\text{F}}$ can be adequately taken as
$0.95$ and $1$, respectively[34, 36]. Then, $\beta$ can be extracted from Eq.
19 after plotting experimental data in FN coordinates and applying previous
procedure. By substituting the extracted $\beta$ in Eq. 18 and using
approximation of the Nordheim function given by Eq. 20, $j_{\text{MG}}$ in Eq.
17 can be calculated. The corresponding theoretical current-density curves are
shown by cyan dots in Fig. 10. It is obvious that both elementary FN equation
(Eq. 15) and MG equation (Eq. 17) fail to explain the saturation trend
observed experimentally.

The discrepancy between theory and experiment is due to the fact that FN type
equations are derived for the metal where the number of carriers is
approaching infinity, surface is equipotential, and resistance on the current
path can be ignored. However, this assumption is obviously not valid for non-
metal and semi-metal emitters: compare the carrier concentration $\sim
10^{23}$ cm-3 in copper to the carrier concentration $\sim 10^{18}$-$10^{19}$
cm-3 in CNT and UNCD. The development of improved theoretical models that
would take semiconductor effects into account is critical. Nevertheless, most
studies relate this saturation effect to impurities at emitter apex[13],
space-charge effect[14, 37], series contact resistance between substrate and
the emitter[38, 39], etc. As we stated before, the space charge is not a
dominant effect at the current-density levels of our interest. Surface
adsorbates do not cause saturation for metal emitters, so we believe it is not
an important effect for our samples. Measurements[16] show that voltage drop
due to series contact resistance is small for typical range of applied
potential, so it can also be ruled out. On the other hand, it is proven both
experimentally[16] and analytically[40] that small voltage drop along the
emitter body can cause large drop in enhancement factor at the emitter apex.
The voltage loss along the body can also be modeled by penetration of the
field into the sample through defining a certain depletion length of $W$ from
the emitter surface[41] into the bulk. Such a $W$ is a result of poor
screening due to limited number of charge carriers in non-metallic field
emitters. Therefore, the effects of the decrease in field enhancement factor
as a result of voltage loss along the emitter[40, 16] or limited charge
supply[41] appear to be better hypotheses. The saturation current density,
i.e. the maximum current density that can be extracted from a sample assuming
carriers reach saturation velocity, is given by[41]

$j_{\text{sat}}=\frac{en^{2/3}v_{\text{sat}}}{W}=\frac{en^{2/3}}{\tau},$ (23)

where $j_{\text{sat}}~{}(\frac{\text{A}}{\text{cm}^{2}})$ is the saturation
current density, $e=1.6\times 10^{-19}~{}\text{C}$ is the electron charge,
$n~{}(\text{cm}^{-3})$ is the electron carrier density,
$v_{\text{sat}}~{}(\frac{\text{cm}}{\text{s}})$ is the saturation velocity,
$W~{}(\text{cm})$ is the length of the sub-surface depletion region, and
$\tau$ is the transit time, i.e. the time it takes for charges to cross over
the depletion region and reach the cathode surface to get emitted into vacuum.
Typical parameters for CNT fiber[2] are: $n\sim 10^{18}~{}\text{cm}^{-3}$,
$v_{\text{sat}}\sim 10^{7}~{}\text{cm}/\text{s}$, and $W\sim 8.9\times
10^{-5}~{}\text{cm}~{}(890~{}\text{nm})$. Then, Eq. 23 gives $1.8\times
10^{4}~{}\frac{\text{A}}{\text{cm}^{2}}$ which is shown with orange dashed
line in Figs. 10B and 10C for reference. Another dashed orange line $3.7\times
10^{5}~{}\frac{\text{A}}{\text{cm}^{2}}$ is also drawn for UNCD in Fig. 10A,
where $n$ is taken as $10^{19}~{}\text{cm}^{-3}$ and $W$ is replaced with the
film thickness of 200 nm atop of a metal substrate that would terminate
electric field lines. Yet another, orange solid line, is presented in the same
plot. This additional line is based on the extended Stratton-Baskin-Lvov-
Fursey (SBLF) theory formulated in our previous work (see Ref. 26). In SBLF
theory, Poisson’s and Stratton’s equations are solved self-consistently at the
material-vacuum interface. The extended SBLF takes into account detailed
density of states of a material and, importantly, accounts for the effect of
field-induced reduction of electron mobility (referred in this work as the
transit time limitation) that plays a central role. The comparison between
solid and dashed orange lines suggests that Eq. 23 can be safely used in place
of detailed and time consuming SBLF calculations if a fast benchmarking of
experimental data is needed.

Fig. 10 demonstrates the behaviour of real field emitters that is
significantly different from that dictated by the classical field emission
theory. First of all, we find that neither Fowler-Nordheim nor Murphy-Good
formulation of $j$-$E$ relation captures the functional behaviour of the
studied field emitters, while demonstrating significant disagreement with
experimental measurements. Thus, the classical field emission appears to have
no relevance toward describing experimental behaviour of this class of
cathodes (be it a first-time turn-on or a consecutive conditioned operation).
Additionally, the 1D vacuum space-charge effect (it is obtained from a 1D
Child-Langmuir law treatment[27] and thus provides the lowest boundary for the
space-charge saturation limit as compared to a 2D case[42]) still
overestimates the experimentally observed saturation by 1 to 3 orders of
magnitude.

The transit-time limited field emission model appears to capture the
experimentally observed emission characteristics unique to the advanced
semiconducting field emission materials, UNCD and CNT. It also correctly
explains two major facts important for cathode conditioning (’Fiber1’ and
’Fiber 2’ case studies in this work). First, despite the characteristic field
emission parameters (e.g., field enhancement factor and turn-on field) may
drastically alter, the current density saturates due to intrinsic material
properties and thus remains the same as long as the material itself does not
alter during conditioning. Second, it correctly implies that once a
semiconductor is depleted, the increase of the output current is only possible
if the emission area increases (see the results in Fig. 9). The emission area
grows with the electric field and so does the output current, but the current
density is a constant of the electric field. This creates a sizeable
misconception when measured $I$-$V$ or $I$-$E$ curves, and not $j$-$E$ curves,
are analyzed using the elementary Fowler-Nordheim approach.

##  IV Conclusions

In conclusion, we demonstrated that the experimental field emission current
density of UNCD and CNT materials remained constant with applied electric
field demonstrating saturation behaviour. These results emphasized that
classical Fowler-Nordheim and Murphy-Good theories fail to describe
fundamental field emission properties of these materials. Moreover, it was
shown that the saturation level cannot be explained by such external effects
like space charge, surface adsorbates, or contact resistance. Instead, it was
found that the transit-time limited field emission, which accounts for the
limited carrier concentration and the saturation drift velocity, provides the
current-density saturation levels very close to the upper limits of
experimentally measured ones for all cases studied in this work.

The direct benchmarking between theory and experiment became possible solely
thanks to the development of a new pattern recognition algorithm. Application
of this algorithm was demonstrated and emphasized by counting the number of
field emission sites and calculating the apparent field emission area from dc
micrographs. The algorithm is fast. Run on a personal laptop with Intel i5 2.5
GHz core, the full set of calculations took 23.8 seconds for Fig. 1, 18.9
seconds for Fig. 4, and 26.3 seconds for Fig. 6. It took the longest
processing time for Fig. 6 because a soft decision boundary was used and false
LMs had to be sorted out.

Future work is being focused on unsupervised machine learning, especially one-
class support vector machines, to replace the Gaussian decision boundary and
thus completely automate the pattern recognition. The current/future releases
are/will be available as a freeware on GitHub[11].

## Acknowledgments

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of High Energy Physics under Award No. DE-SC0020429.
O.C. was supported by the U.S. National Science Foundation under Award
PHY-1549132, the Center for Bright Beams. We would like to thank Prof. Richard
Forbes for his constructive critique and valuable suggestions.

## References

  * Chubenko  _et al._ [2017] O. Chubenko, S. S. Baturin, K. K. Kovi, A. V. Sumant, and S. V. Baryshev, ACS Applied Materials & Interfaces 9, 33229 (2017).
  * Posos  _et al._ [2020] T. Y. Posos, S. B. Fairchild, J. Park, and S. V. Baryshev, Journal of Vacuum Science & Technology B 38, 024006 (2020).
  * Shao  _et al._ [2016] J. Shao, J. Shi, S. P. Antipov, S. V. Baryshev, H. Chen, M. Conde, W. Gai, G. Ha, C. Jing, F. Wang, and E. Wisniewski, Phys. Rev. Lett. 117, 084801 (2016).
  * Shao  _et al._ [2019] J. Shao, M. Schneider, G. Chen, T. Nikhar, K. K. Kovi, L. Spentzouris, E. Wisniewski, J. Power, M. Conde, W. Liu, and S. V. Baryshev, Phys. Rev. Accel. Beams 22, 123402 (2019).
  * Cui and Robertson [2002] J. B. Cui and J. Robertson, Journal of Vacuum Science & Technology B 20, 238 (2002).
  * Popov  _et al._ [2018] E. O. Popov, A. G. Kolosko, S. V. Filippov, and E. I. Terukov, Journal of Vacuum Science & Technology B 36, 02C106 (2018).
  * Serbun  _et al._ [2013] P. Serbun, B. Bornmann, A. Navitski, G. Müller, C. Prommesberger, C. Langer, F. Dams, and R. Schreiner, Journal of Vacuum Science & Technology B 31, 02B101 (2013).
  * Schreiner  _et al._ [2013] R. Schreiner, C. Langer, C. Prommesberger, S. Mingels, P. Serbun, and G. Müller, in  _2013 26th International Vacuum Nanoelectronics Conference (IVNC)_ (2013) pp. 1–2.
  * Köck  _et al._ [2004] F. Köck, J. Garguilo, and R. Nemanich, Diamond and Related Materials 13, 1022 (2004).
  * Lyashenko  _et al._ [2009] S. A. Lyashenko, A. P. Volkov, R. R. Ismagilov, and A. N. Obraztsov, Technical Physics Letters 35, 249 – 252 (2009).
  * Git [2021] FEPIC, https://github.com/schne525/FEPic (2021).
  * Perry [1962] R. L. Perry, J. Appl. Phys. 33, 1875 (1962).
  * Dean and Chalamala [2000] K. A. Dean and B. R. Chalamala, Applied Physics Letters 76, 375 (2000).
  * Cahay  _et al._ [2014] M. Cahay, P. T. Murray, T. C. Back, S. Fairchild, J. Boeckl, J. Bulmer, K. K. K. Koziol, G. Gruen, M. Sparkes, F. Orozco, and W. O’Neill, Applied Physics Letters 105, 173107 (2014).
  * Serbun [2014] P. Serbun,  _A systematic investigation of carbon, metallic and semiconductor nanostructures for field-emission cathode applications_ , Ph.D. thesis, University of Wuppertal (2014).
  * Minoux  _et al._ [2005] E. Minoux, O. Groening, K. B. K. Teo, S. H. Dalal, L. Gangloff, J.-P. Schnell, L. Hudanski, I. Y. Y. Bu, P. Vincent, P. Legagneux, G. A. J. Amaratunga, and W. I. Milne, Nano Letters 5, 2135 (2005).
  * Ducati  _et al._ [2002] C. Ducati, E. Barborini, P. Piseri, P. Milani, and J. Robertson, J. Appl. Phys. 92, 5482 (2002).
  * Rodriguez and Laio [2014] A. Rodriguez and A. Laio, Science 344, 1492 (2014).
  * Verbeek  _et al._ [1988] P. Verbeek, H. Vrooman, and L. Van Vliet, Signal Processing 15, 249 (1988).
  * Mukherjee  _et al._ [2020] D. Mukherjee, L. Miao, G. Stone, and N. Alem, Advanced Structural and Chemical Imaging 6, 1 (2020).
  * Lindeberg [1998] T. Lindeberg, International Journal of Computer Vision 30, 79 (1998).
  * Fowler and Nordheim [1928] R. H. Fowler and L. Nordheim, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 119, 173 (1928).
  * Chen  _et al._ [2019] G. Chen, G. Adhikari, L. Spentzouris, K. K. Kovi, S. Antipov, C. Jing, W. Andreas Schroeder, and S. V. Baryshev, Applied Physics Letters 114, 093103 (2019).
  * Nikhar  _et al._ [2020] T. Nikhar, S. V. Baryshev, G. Adhikari, and A. W. Schroeder, arXiv (2020), arXiv:2011.00722 [physics.acc-ph] .
  * Filip  _et al._ [2001] V. Filip, D. Nicolaescu, M. Tanemura, and F. Okuyama, Ultramicroscopy 89, 39 (2001), 46th International Field Emission Symposium.
  * Chubenko  _et al._ [2019] O. Chubenko, S. S. Baturin, and S. V. Baryshev, Journal of Applied Physics 125, 205303 (2019).
  * Barbour  _et al._ [1953] J. P. Barbour, W. W. Dolan, J. K. Trolan, E. E. Martin, and W. P. Dyke, Phys. Rev. 92, 45 (1953).
  * Sommerfeld and Bethe [1934] H. Sommerfeld and H. Bethe, Handbuch der Physik 24, 436 (1934).
  * Murphy and Good [1956] E. L. Murphy and R. H. Good, Phys. Rev. 102, 1464 (1956).
  * Houston [1952] J. M. Houston, Phys. Rev. 88, 349 (1952).
  * Forbes [2019] R. G. Forbes, Journal of Applied Physics 126, 210901 (2019).
  * Nordhiem and Fowler [1928] L. W. Nordhiem and R. H. Fowler, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 121, 626 (1928).
  * Fursey [2005] G. N. Fursey,  _Field Emission in Vacuum Microelectronics_, Microdevices (Springer US, 2005).
  * Burgess  _et al._ [1953] R. E. Burgess, H. Kroemer, and J. M. Houston, Phys. Rev. 90, 515 (1953).
  * Forbes [2006] R. G. Forbes, Applied Physics Letters 89, 113122 (2006).
  * Forbes [2020] R. G. Forbes, Review of Scientific Instruments 91, 107101 (2020).
  * Feng and Verboncoeur [2006] Y. Feng and J. P. Verboncoeur, Physics of Plasmas 13, 073105 (2006).
  * Luginsland  _et al._ [1996a] J. W. Luginsland, A. Valfells, and Y. Y. Lau, Applied Physics Letters 69, 2770 (1996a).
  * Zhang  _et al._ [2017] P. Zhang, S. B. Fairchild, T. C. Back, and Y. Luo, AIP Advances 7, 125203 (2017).
  * Forbes [2017] R. G. Forbes, Applied Physics Letters 110, 133109 (2017).
  * Baturin  _et al._ [2018] S. S. Baturin, A. V. Zinovev, and S. V. Baryshev, Current saturation in nonmetallic field emitters (2018), arXiv:1710.03692 [cond-mat.mes-hall] .
  * Luginsland  _et al._ [1996b] J. W. Luginsland, Y. Y. Lau, and R. M. Gilgenbach, Phys. Rev. Lett. 77, 4668 (1996b).