# Distributed Gauss-Newton Method for AC State Estimation: A Belief
Propagation Approach

Mirsad Cosovic  Schneider Electric DMS NS LLC
Novi Sad, Serbia
Email<EMAIL_ADDRESS>Dejan Vukobratovic
Department of Power, Electronics and Communications Engineering,
University of Novi Sad, Serbia
Email<EMAIL_ADDRESS>
###### Abstract

In this paper, we propose a solution to an AC state estimation problem in
electric power systems using a fully distributed Gauss-Newton method. The
proposed method is placed within the context of factor graphs and belief
propagation algorithms and closed-form expressions for belief propagation
messages exchanged along the factor graph are derived. The obtained algorithm
provides the same solution as the conventional weighted least-squares state
estimation. Using a simple example, we provide a step-by-step presentation of
the proposed algorithm. Finally, we discuss the convergence behaviour using
the IEEE 14 bus test case.

###### Index Terms:

AC State Estimation, Electric Power System, Factor Graphs, Gaussian Belief
Propagation, Distributed Gauss-Newton Method

##  I Introduction

The state estimation (SE) is an essential part of the real-time energy
management system (EMS), and it provides inputs for other EMS functions. A SE
algorithm, jointly with network topology processors, observability analysis
and bad data analysis, provides an estimate of the system state according to
the network topology and available measurements. SE is performed on a
bus/branch model and used to reconstruct the state of the observable part of
the system. Conventional SE algorithms are centralized and typically use the
Gauss-Newton method to solve the non-linear weighted least-squares (WLS)
problem [1], [2].

Deregulation of electric power systems implies their decentralized structure,
however, integrated control and monitoring across the entire network is still
needed. In view of recent trends, control centers are mostly migrating toward
distributed control centers [3]. Consequently, many centralized algorithms of
EMS have to be redefined, requirements being distributed and computationally
more efficient algorithms.

Probabilistic graphical models, such as factor graphs (FGs), represent a
powerful tool for modeling probabilistic systems. The algorithm for exact
inference on probabilistic graphical models without loops is known as the
belief propagation (BP) algorithm [4], [5]. Using BP, it is possible to
efficiently calculate marginal distributions or a mode of the joint
distribution of the system of random variables. The BP algorithm can be also
applied to graphical models with loops (loopy BP)[6], although in that case,
the solution is not guaranteed to converge to the correct marginals/modes of
the joint distribution. BP is a fully distributed algorithm that takes
probability distributions as an input, processes them, and outputs marginal
probability distributions used to estimate values of state variables. This
makes it a flexible solution for accommodation of distributed power sources
and time-varying loads in various applications of electric power systems.

The work in [7] provides the first demonstration of BP applied to the SE
problem. Although this work is elaborate in terms of using, e.g.,
environmental correlation via historical data, it applies BP to a simple
linearized DC model. The AC model is recently addressed in [8], where tree-
reweighted BP is applied using preprocessed weights obtained by randomly
sampling the space of spanning trees.

In this paper, we also solve the AC model via BP but in a completely different
framework that we find simpler and more intuitive. We consider the AC SE model
that we cast into a FG representation and solve using the BP algorithm. The
proposed BP algorithm is obtained after a transformation of the initial WLS
problem into the maximum a posteriori probability (MAP) problem that estimates
the vector of increments of the state variables. The resulting BP algorithm
has the interpretation of a fully distributed Gauss-Newton method with the
same accuracy as the conventional or centralized SE. Consequently, the
presented algorithm applies similar mathematical framework as the conventional
SE, which simplifies its integration into existing EMS. In addition, it can be
implemented in a fully distributed manner suitable for the distributed multi-
area SE environment.

The structure of this paper is as follows: Section II describes the
conventional (centralized) SE and defines an optimization problem which allows
a solution via the BP approach. Section III formulates closed form expressions
for BP messages. In Section IV, we give a step-by-step description of the
proposed algorithm, while Section V considers the convergence performance and
numerical results for the IEEE 14 bus test case. Concluding remarks are
included in Section VI.

##  II Electric power system state estimation

The AC SE problem reduces to solving the system of equations[9]:

$\displaystyle\mathbf{z}=\mathbf{h}(\mathbf{x})+\mathbf{u},$ (1)

where $\mathbf{h}(\mathbf{x})=(h_{1}(\mathbf{x}),\dots,h_{k}(\mathbf{x}))$
includes both non-linear and linear measurement functions (see Appendix for
details), $\mathbf{x}=(x_{1},\dots,x_{n})$ is the vector of the state
variables, $\mathbf{z}=(z_{1},\dots,z_{k})$ is the vector of independent
measurements (where $n<k$), and $\mathbf{u}=(u_{1},\dots,u_{k})$ is the vector
of measurement errors.

The state variables are bus voltage magnitudes and bus voltage angles,
transformer magnitudes of turns ratio and transformer angles of turns ratio.
Without loss of generality, in the rest of the paper, we observe bus voltage
angles $\bm{\uptheta}=$ $(\theta_{1},\dots,\theta_{N})$ and bus voltage
magnitudes $\mathbf{V}=$ $(V_{1},\dots,V_{N})$ as state variables
$\mathbf{x}\equiv(\bm{\uptheta},\mathbf{V})$, where $N$ is the number of buses
($n=2N$).

Under the assumption that measurement errors $\mathbf{u}$ follow a zero-mean
Gaussian distribution, the probability density function associated with the
m-th measurement equals:

$\begin{gathered}\mathcal{N}(z_{m}|\mathbf{x},\sigma_{m}^{2})=\cfrac{1}{\sigma_{m}\sqrt{2\pi}}\exp\Bigg{\\{}\cfrac{[z_{m}-h_{m}(\mathbf{x})]^{2}}{2\sigma_{m}^{2}}\Bigg{\\}},\end{gathered}$
(2)

where $z_{m}$ is the value of the measurement, $\sigma_{m}^{2}$ is the
measurement variance, and the function $h_{m}(\mathbf{x})$ connects the vector
of state variables to the value of m-th measurement.

One can find the MAP solution to the SE problem via maximization of the
likelihood function, which is defined via likelihoods of $k$ independent
measurements:

$\begin{gathered}\hat{\mathbf{x}}=\mathrm{arg}\max_{\mathbf{x}}\mathcal{L}(\mathbf{z}|\mathbf{x})=\mathrm{arg}\max_{\mathbf{x}}\prod_{h=1}^{k}\mathcal{N}(z_{h}|\mathbf{x},\sigma_{h}^{2}).\end{gathered}$
(3)

The maximum likelihood estimator (3) is equivalent to the weighted least-
squares estimator whose solution can be found using the Gauss-Newton method:

$\displaystyle\mathbf{J}(\mathbf{x}^{\nu})^{\mathrm{T}}\mathbf{W}\mathbf{J}(\mathbf{x}^{\nu})\Delta\mathbf{x}^{\nu}=\mathbf{J}(\mathbf{x}^{\nu})^{\mathrm{T}}\mathbf{W}\mathbf{r}(\mathbf{x}^{\nu})$
(.1)
$\displaystyle\mathbf{x}^{\nu+1}=\mathbf{x}^{\nu}+\Delta\mathbf{x}^{\nu},$
(.2)

where $\nu$ is the iteration step, $\Delta\mathbf{x}^{\nu}\in\mathbb{R}^{n}$
is the vector of increments of the state variables,
$\mathbf{J}(\mathbf{x}^{\nu})\in\mathbb{R}^{k\mathrm{x}n}$ is the Jacobian
matrix of measurement functions $\mathbf{h}(\mathbf{x}^{\nu})$ (see Appendix
for details), $\mathbf{W}\in\mathbb{R}^{k\mathrm{x}k}$ is a diagonal matrix
containing inverses of measurement variances, and
$\mathbf{r}(\mathbf{x}^{\nu})=\mathbf{z}-\mathbf{h}(\mathbf{x}^{\nu})$ is the
vector of residuals[1].

At each iteration $\nu$, the Gauss-Newton method returns a new estimate of
$\mathbf{x}$, which in a given iteration may be observed as a constant vector.
If the Jacobian matrix $\mathbf{J}(\mathbf{x}^{\nu})$ has a full column rank,
the equation (.1) represents the linear WLS solution of the minimization
problem [10]:

$\begin{gathered}\min_{\Delta\mathbf{x}^{\nu}}||\mathbf{W}^{1/2}[\mathbf{r}(\mathbf{x}^{\nu})-\mathbf{J}(\mathbf{x}^{\nu})\Delta\mathbf{x}^{\nu}]||_{2}^{2}.\end{gathered}$
(5)

Hence, at each iteration $\nu$, we can consider system of linear equations:

$\displaystyle\mathbf{r}(\mathbf{x}^{\nu})=\mathbf{g}(\Delta\mathbf{x}^{\nu})+\mathbf{u},$
(6)

where
$\mathbf{g}(\Delta\mathbf{x}^{\nu})=\mathbf{J}(\mathbf{x}^{\nu})\Delta\mathbf{x}^{\nu}$
comprises linear functions. The equation (.1) is the weighted normal equation
for the minimization problem defined as (5), or alternatively, equation (.1)
is WLS solution of (6).

Consequently, the probability density function associated with the m-th
measurement (i.e., the m-th residual component $r_{m}$) at any iteration
$\nu$:

$\begin{gathered}\mathcal{N}(r_{m}(\mathbf{x}^{\nu})|{\Delta\mathbf{x}^{\nu}},\sigma_{m}^{2})=\cfrac{1}{\sigma_{m}\sqrt{2\pi}}\exp\Bigg{\\{}\cfrac{[r_{m}(\mathbf{x}^{\nu})-g_{m}(\Delta\mathbf{x}^{\nu})]^{2}}{2\sigma_{m}^{2}}\Bigg{\\}}.\end{gathered}$
(7)

The MAP solution of (3) can be redefined as an iterative optimization problem
where, instead of solving (.1) and (.2), we solve MAP (sub)problem:

$\displaystyle\Delta\hat{\mathbf{x}}^{\nu}$
$\displaystyle=\mathrm{arg}\max_{\Delta\mathbf{x}^{\nu}}\mathcal{L}\Big{(}\mathbf{r}(\mathbf{x}^{\nu})|\Delta\mathbf{x}^{\nu}\Big{)}$
(8)
$\displaystyle=\mathrm{arg}\max_{\Delta\mathbf{x}^{\nu}}\prod_{h=1}^{k}\mathcal{N}\Big{(}r_{h}(\mathbf{x}^{\nu})|\Delta\mathbf{x}^{\nu},\sigma_{h}^{2}\Big{)}$
$\displaystyle\mathbf{x}^{\nu+1}$
$\displaystyle=\mathbf{x}^{\nu}+\Delta\hat{\mathbf{x}}^{\nu}.$

As we show next, the solution to the above MAP subproblem over increment
variables $\Delta{\mathbf{x}}^{\nu}$ can be efficiently obtained using BP
algorithm applied over the underlying factor graph.

Note that, if the factor graph corresponding to the problem (8) (see Section
III) is a tree, the resulting BP algorithm provides a solution equal to the
linear WLS solution $\Delta{\mathbf{x}}^{\nu}$ of (.1). In general, if the
factor graph contains loops, the BP solution of $\Delta\hat{\mathbf{x}}^{\nu}$
in each iteration $\nu$ (outer iteration loop) will be obtained via iterative
BP algorithm (inner iteration loops). Every inner BP iteration loop
$\rho=1,2,\dots,\tau(\nu)$ outputs $\Delta\hat{\mathbf{x}}^{\nu,\rho}$, where
$\tau(\nu)$ is the number of inner BP iterations within outer iteration $\nu$.

##  III Factor Graphs and BP algorithm

Factor graphs and BP algorithm are widely used tools for probabilistic
inference [4], [5]. In our scenario, FGs consist of variable nodes for every
variable in the likelihood function and of factor nodes for each likelihood
factor. Therefore, for the MAP subproblem defined in (8), the increments
$\Delta\mathbf{x}$ of state variables $\mathbf{x}$ will appear as variable
nodes, while residuals $\mathbf{r}$ of measurements $\mathbf{z}$ will define
factor nodes. A factor node connects to a variable node if the increment
variable is an argument of the corresponding function
${g_{m}}({\Delta\mathbf{x}})$, which is equivalent to say that the
corresponding state variable is an argument of the measurement function
$h_{m}(\mathbf{x})$.

Consider the part of a factor graph shown in Fig. 1 with group of factor nodes
$\mathcal{F}=\\{f_{s},f_{w},...,f_{W}\\}$ that are neighbours of the variable
node $\Delta x_{m}$.

Figure 1: Message from variable node $\Delta x_{m}$ to factor node $f_{s}$

The message from the variable node $\Delta x_{m}$ to the factor node $f_{s}$
is equal to the product of all incoming factor node to variable node messages
arriving at all the other incident edges:

$\begin{gathered}\mu_{\Delta x_{m}\to f_{s}}(\Delta
x_{m})=\prod_{f_{a}\in\mathcal{F}\setminus f_{s}}\mu_{f_{a}\to\Delta
x_{m}}(\Delta x_{m}),\end{gathered}$ (9)

where $\mathcal{F}\setminus f_{s}$ defines the set of factor nodes incident to
the variable node $\Delta x_{m}$, excluding the factor node $f_{s}$.

It can be shown that the message $\mu_{\Delta x_{m}\to f_{s}}(\Delta x_{m})$
is represented by the Gaussian function:

$\begin{gathered}\mu_{\Delta x_{m}\to f_{s}}(\Delta
x_{m})\propto\mathcal{N}(r_{\Delta x_{m}\to f_{s}}|\Delta x_{m},\sigma_{\Delta
x_{m}\to f_{s}}^{2}),\end{gathered}$ (10)

with mean $r_{\Delta x_{m}\to f_{s}}$ and variance $\sigma_{\Delta x_{m}\to
f_{s}}^{2}$:

$\displaystyle r_{\Delta x_{m}\to f_{s}}$
$\displaystyle=\Bigg{(}\sum_{f_{a}\in\mathcal{F}\setminus
f_{s}}\cfrac{r_{f_{a}\to\Delta x_{m}}}{\sigma_{f_{a}\to\Delta
x_{m}}^{2}}\Bigg{)}\sigma_{\Delta x_{m}\to f_{s}}^{2}$ (11)
$\displaystyle\cfrac{1}{\sigma_{\Delta x_{m}\to f_{s}}^{2}}$
$\displaystyle=\sum_{f_{a}\in\mathcal{F}\setminus
f_{s}}\cfrac{1}{\sigma_{f_{a}\to\Delta x_{m}}^{2}}.$

Consider the part of the factor graph shown in Fig. 2 that consists of the
group of variable nodes $\mathcal{X}\in\\{\Delta x_{m},\Delta x_{l},...,\Delta
x_{L}\\}$ that are neighbours of the factor node $f_{s}$. The message from the
factor node $f_{s}$ to the variable node $\Delta x_{m}$ is defined as a
product of all incoming variable node to factor node messages arriving at all
the other incident edges multiplied by the Gaussian function associated to the
factor node $f_{s}$ and marginalized over all of the variables associated with
the incoming messages:

$\displaystyle\mu_{f_{s}\to\Delta x_{m}}(\Delta
x_{m})=\int\displaylimits_{\Delta x_{l}}\dots\int\displaylimits_{\Delta
x_{L}}\mathcal{N}(r_{f_{s}}|\Delta x_{m},\Delta x_{l}\dots\Delta
x_{L},\sigma_{f_{s}}^{2})$ (12) $\displaystyle\prod_{\Delta
x_{b}\in\mathcal{X}\setminus\Delta x_{m}}\mu_{\Delta x_{b}\to f_{s}}(\Delta
x_{b})\cdot\mathrm{d}\Delta x_{b},$

where $\mathcal{X}\setminus\Delta x_{m}$ is the set of variable nodes incident
to the factor node $f_{s}$, excluding the variable node $\Delta x_{m}$.

Figure 2: Message from factor node $f_{s}$ to variable node $\Delta x_{m}$

It can be shown that the message $\mu_{f_{s}\to\Delta x_{m}}(\Delta x_{m})$ is
represented by the Gaussian function:

$\begin{gathered}\mu_{f_{s}\to\Delta x_{m}}(\Delta
x_{m})\propto\mathcal{N}(r_{f_{s}\to\Delta x_{m}}|\Delta
x_{m},\sigma_{f_{s}\to\Delta x_{m}}^{2}),\end{gathered}$ (13)

with mean $r_{f_{s}\to\Delta x_{m}}$ and variance $\sigma_{f_{s}\to\Delta
x_{m}}^{2}$:

$\begin{gathered}r_{f_{s}\to\Delta x_{m}}=\cfrac{1}{C_{\Delta
x_{m}}}\Bigg{(}r_{f_{s}}-\sum_{\Delta x_{b}\in\mathcal{X}\setminus\Delta
x_{m}}C_{\Delta x_{b}}\cdot r_{\Delta x_{b}\to f_{s}}\Bigg{)}\\\
\sigma_{f_{s}\to\Delta x_{m}}^{2}=\cfrac{1}{C_{\Delta
x_{m}}^{2}}\Bigg{(}\sigma_{f_{s}}^{2}+\sum_{\Delta
x_{b}\in\mathcal{X}\setminus\Delta x_{m}}C_{\Delta
x_{b}}^{2}\cdot\sigma_{\Delta x_{b}\to f_{s}}^{2}\Bigg{)}.\end{gathered}$ (14)

The coefficients $C_{\Delta x_{i}},\;i=m,l\dots,L$, are Jacobian elements of
the measurement function (see Appendix for details) associated with the factor
node $f_{s}$:

$\begin{gathered}C_{\Delta x_{i}}=\cfrac{\partial
h(x_{m},x_{l},\dots,x_{L})}{\partial x_{i}},\;i=m,l\dots,L.\end{gathered}$
(15)

Note that, due to the fact that all the measurements follow Gaussian
distribution and that BP processing in both variable and function nodes
preserve ”Gaussianity”, all the messages exchanged in the presented BP are
Gaussian distributions. The resulting BP algorithm is known as Gaussian BP
algorithm in which all the BP messages are completely represented using only
means and variances[11].

The marginal of state variable increment $\Delta x_{m}$, illustrated in Fig.
3, is obtained as the product of all incoming messages into the variable node
$\Delta x_{m}$:

$\begin{gathered}p(\Delta
x_{m})=\prod_{f_{c}\in\mathcal{F}}\mu_{f_{c}\to\Delta x_{m}}(\Delta
x_{m}),\end{gathered}$ (16)

where $\mathcal{F}$ is the set of factor nodes incident to the variable node
$\Delta x_{m}$.

Figure 3: Marginal inference for $\Delta x_{m}$

Thus the marginal has Gaussian form:

$\begin{gathered}p(\Delta x_{m})\propto\mathcal{N}(\Delta\hat{x}_{m}|\Delta
x_{m},\sigma_{\Delta x_{m}}^{2}),\end{gathered}$ (17)

with mean $\Delta\hat{x}_{m}$ which represents the estimated value of the
state variable increment $\Delta x_{m}$ and variance $\sigma_{\Delta
x_{m}}^{2}$:

$\displaystyle\Delta\hat{x}_{m}$
$\displaystyle=\Bigg{(}\sum_{f_{c}\in\mathcal{F}}\cfrac{r_{f_{c}\to\Delta
x_{m}}}{\sigma_{f_{c}\to\Delta x_{m}}^{2}}\Bigg{)}\sigma_{\Delta x_{m}}^{2}$
(18) $\displaystyle\cfrac{1}{\sigma_{\Delta x_{m}}^{2}}$
$\displaystyle=\sum_{f_{c}\in\mathcal{F}}\cfrac{1}{\sigma_{f_{c}\to\Delta
x_{m}}^{2}}.$

The MAP subproblem defined in (8) can be efficiently solved using (11), (14)
and (18).

##  IV Toy Example

An illustrative example presented in Fig. 4 will be used to provide a step-by-
step presentation of the proposed algorithm. The simple three bus radial
network contains three direct measurement devices that directly measure state
variables $M_{dir}\in\\{M_{V_{1}},M_{\theta_{2}},M_{\theta_{3}}\\}$, and two
indirect measurement devices $M_{ind}\in\\{M_{P_{12}},M_{P_{23}}\\}$ that
measure state variables indirectly.

Figure 4: Bus/branch model

Input data for SE from measurement devices are Gaussian-type functions
represented by means and variances:
$z_{dir}\in\\{z_{V_{1}},z_{\theta_{2}},z_{\theta_{3}}\\}$,
$\sigma_{dir}^{2}\in\\{\sigma_{V_{1}}^{2},\sigma_{\theta_{2}}^{2},\sigma_{\theta_{3}}^{2}\\}$
and $z_{ind}\in\\{z_{P_{12}},z_{P_{23}}\\}$ ,
$\sigma_{ind}^{2}\in\\{\sigma_{P_{12}}^{2},\sigma_{P_{23}}^{2}\\}$.

The corresponding FG is given in Fig. 5, where we define indirect factor nodes
$f_{r_{P_{12}}}$, $f_{r_{P_{23}}}$ (orange squares) corresponding to indirect
measurements and four types of singly-connected factor nodes (local factor
nodes) described next. The slack factor node $f_{r_{\theta_{1}}}$ (yellow
square) corresponds to the slack or reference bus where the voltage angle has
a given value, therefore, the residual of that state variable is equal to
zero. The direct factor nodes $f_{r_{V_{1}}}$, $f_{r_{\theta_{2}}}$,
$f_{r_{\theta_{3}}}$ (red squares) correspond to the direct measurements. The
initialization factor node $f_{r_{V_{2}}}$ (green square) is needed to start
the algorithm, while the virtual factor node $f_{r_{V_{3}}}$ (blue square) is
used to form a message from a variable node to a factor node. In general, if
the variable node is not directly measured and is singly-connected to the rest
of the FG, we attach a virtual factor node to this variable node.

Figure 5: Factor graph of the illustrative example

In the following, for notational convenience, we denote the variance as
follows: $\sigma^{2}\equiv\xi$.

Algorithm Initialization

  1. 1.

The AC SE in electric power systems assumes ”flat start” or a priori given
values of state variables:

$\begin{gathered}\mathbf{x}^{\\{\nu=0\\}}=[\theta_{1}\;\theta_{2}\;\theta_{3}\;V_{1}\;V_{2}\;V_{3}]^{\\{\nu=0\\}}.\end{gathered}$

  2. 2.

The value of the slack factor node is set to $r_{\theta_{1}}=0$ with variance
$\xi_{\theta_{1}}\to 0$.

  3. 3.

The value of initialization factor nodes and virtual factor nodes are set to
$r_{V_{2}}\to 0$ and $r_{V_{3}}\to 0$, with variances $\xi_{V_{2}}\to\infty$
and $\xi_{V_{3}}\to\infty$.

Iterate - Outer Loop: $\bm{\nu}\mathbf{=0,1,2,\dots};\bm{\rho}\mathbf{=0}$

  1. 4.

Each direct factor node computes residual, e.g.:

$\begin{gathered}r_{\theta_{2}}^{\\{\nu\\}}=z_{\theta_{2}}-\theta_{2}^{\\{\nu\\}}.\end{gathered}$

  2. 5.

Local factor nodes send messages represented by a triplet (residual, variance,
state variable), to incident variable nodes, e.g.:

$\begin{gathered}\mu_{f_{r_{\theta_{1}}}\to\Delta\theta_{1}}^{\\{\nu\\}}:=\big{(}r_{\theta_{1}},\xi_{\theta_{1}},\theta_{1}^{\\{\nu\\}}\big{)}\\\
\mu_{f_{r_{\theta_{2}}}\to\Delta\theta_{2}}^{\\{\nu\\}}:=\big{(}r_{\theta_{2}}^{\\{\nu\\}},\xi_{\theta_{2}},\theta_{2}^{\\{\nu\\}}\big{)}.\end{gathered}$

  3. 6.

Variable nodes forward the incoming messages received from local factor nodes
along remaining edges, e.g.:

$\displaystyle\mu_{\Delta\theta_{2}\to f_{r_{P_{12}}}}^{\\{\nu\\}\\{\rho\\}}$
$\displaystyle=\mu_{f_{r_{\theta_{2}}}\to\Delta\theta_{2}}^{\\{\nu\\}}:=\big{(}r_{\theta_{2}}^{\\{\nu\\}},\xi_{\theta_{2}},\theta_{2}^{\\{\nu\\}}\big{)}$
$\displaystyle:=\big{(}r_{\Delta\theta_{2}\to
f_{r_{P_{12}}}}^{\\{\nu\\}\\{\rho\\}},\xi_{\Delta\theta_{2}\to
f_{r_{P_{12}}}}^{\\{\nu\\}\\{\rho\\}},\theta_{2}^{\\{\nu\\}\\{\rho\\}}\big{)}$
$\displaystyle\mu_{\Delta\theta_{2}\to f_{r_{P_{23}}}}^{\\{\nu\\}\\{\rho\\}}$
$\displaystyle=\mu_{f_{r_{\theta_{2}}}\to\Delta\theta_{2}}^{\\{\nu\\}}:=\big{(}r_{\theta_{2}}^{\\{\nu\\}},\xi_{\theta_{2}},\theta_{2}^{\\{\nu\\}}\big{)}$
$\displaystyle:=\big{(}r_{\Delta\theta_{2}\to
f_{r_{P_{23}}}}^{\\{\nu\\}\\{\rho\\}},\xi_{\Delta\theta_{2}\to
f_{r_{P_{23}}}}^{\\{\nu\\}\\{\rho\\}},\theta_{2}^{\\{\nu\\}\\{\rho\\}}\big{)}.$

  4. 7.

Indirect factor nodes compute residuals, e.g.:

$\begin{gathered}r_{P_{12}}^{\\{\nu\\}}=z_{P_{12}}-h_{P_{12}}(\theta_{1}^{\\{\nu\\}},\theta_{2}^{\\{\nu\\}},V_{1}^{\\{\nu\\}},V_{2}^{\\{\nu\\}}).\end{gathered}$

  5. 8.

Indirect factor nodes compute appropriate Jacobian elements associated with
state variables, e.g.:

$\displaystyle
C_{P_{12},\Delta\theta_{1}}^{\\{\nu\\}}=\cfrac{\mathrm{\partial}{h_{P_{12}}(\cdot)}}{\mathrm{\partial}\theta_{1}}$
$\displaystyle={V}_{1}^{\\{\nu\\}}{V}_{2}^{\\{\nu\\}}(g_{12}\mbox{sin}\theta_{12}^{\\{\nu\\}}-b_{12}\mbox{cos}\theta_{12}^{\\{\nu\\}})$
$\displaystyle C_{P_{12},\Delta
V_{2}}^{\\{\nu\\}}=\cfrac{\mathrm{\partial}{h_{{P_{12}}}(\cdot)}}{\mathrm{\partial}V_{2}}$
$\displaystyle=-{V}_{1}^{\\{\nu\\}}(g_{12}\mbox{cos}\theta_{12}^{\\{\nu\\}}+b_{12}\mbox{sin}\theta_{12}^{\\{\nu\\}}).$

Iterate - Inner Loop: $\bm{\rho}\mathbf{=1,2,\dots,\bm{\tau}(\bm{\nu})}$

  1. 9.

Indirect factor nodes send messages as pairs along incident edges according to
(14), e.g.:

$\displaystyle\mu_{f_{r_{P_{12}}}\to\Delta\theta_{2}}^{\\{\rho\\}}:=\big{(}r_{f_{r_{P_{12}}}\to\Delta\theta_{2}}^{\\{\rho\\}},\xi_{f_{r_{P_{12}}}\to\Delta\theta_{2}}^{\\{\rho\\}}\big{)}$
$\displaystyle
r_{f_{r_{P_{12}}}\to\Delta\theta_{2}}^{\\{\rho\\}}=\cfrac{1}{C_{P_{12},\Delta\theta_{2}}^{\\{\nu\\}}}\Big{[}r_{P_{12}}^{\\{\nu\\}}-C_{P_{12},\Delta\theta_{1}}^{\\{\nu\\}}\cdot
r_{\Delta\theta_{1}\to f_{r_{P_{12}}}}^{\\{\nu\\}\\{\rho-1\\}}$
$\displaystyle-C_{P_{12},\Delta V_{1}}^{\\{\nu\\}}\cdot r_{\Delta V_{1}\to
f_{r_{P_{12}}}}^{\\{\nu\\}\\{\rho-1\\}}C_{P_{12},\Delta
V_{2}}^{\\{\nu\\}}\cdot r_{\Delta V_{2}\to
f_{r_{P_{12}}}}^{\\{\nu\\}\\{\rho-1\\}}\Big{]}$
$\displaystyle\xi_{f_{r_{P_{12}}}\to\Delta\theta_{2}}^{\\{\rho\\}}=\cfrac{1}{(C_{P_{12},\Delta\theta_{2}}^{\\{\nu\\}})^{2}}\Big{[}\xi_{P_{12}}+(C_{P_{12},\Delta\theta_{1}}^{\\{\nu\\}})^{2}\cdot\xi_{\Delta\theta_{1}\to
f_{r_{P_{12}}}}^{\\{\nu\\}\\{\rho-1\\}}$ $\displaystyle+(C_{P_{12},\Delta
V_{1}}^{\\{\nu\\}})^{2}\cdot\xi_{\Delta V_{1}\to
f_{r_{P_{12}}}}^{\\{\nu\\}\\{\rho-1\\}}+({C_{P_{12},\Delta
V_{2}}^{\\{\nu\\}}})^{2}\cdot\xi_{\Delta V_{2}\to
f_{r_{P_{12}}}}^{\\{\nu\\}\\{\rho-1\\}}\Big{]}.$

  2. 10.

Variable nodes send messages as pairs along incident edges to indirect factor
nodes according to (11), e.g.:

$\begin{gathered}\mu_{\Delta\theta_{2}\to
f_{r_{P_{12}}}}^{\\{\nu\\}\\{\rho\\}}:=\big{(}r_{\Delta\theta_{2}\to
f_{r_{P_{12}}}}^{\\{\nu\\}\\{\rho\\}},\xi_{\Delta\theta_{2}\to
f_{r_{P_{12}}}}^{\\{\nu\\}\\{\rho\\}}\big{)}\end{gathered}$
$\displaystyle\cfrac{1}{\xi_{\Delta\theta_{2}\to
f_{r_{P_{12}}}}^{\\{\nu\\}\\{\rho\\}}}$
$\displaystyle=\cfrac{1}{\xi_{\theta_{2}}}+\cfrac{1}{\xi_{f_{r_{P_{23}}}\to\Delta\theta_{2}}^{\\{\rho\\}}}$
$\displaystyle r_{\Delta\theta_{2}\to f_{r_{P_{12}}}}^{\\{\nu\\}\\{\rho\\}}$
$\displaystyle=\Bigg{(}\cfrac{r_{\theta_{2}}^{\\{\nu\\}}}{\xi_{\theta_{2}}}+\cfrac{r_{f_{r_{P_{23}}}\to\Delta\theta_{2}}^{\\{\rho\\}}}{\xi_{f_{r_{P_{23}}}\to\Delta\theta_{2}}^{\\{\rho\\}}}\Bigg{)}\xi_{\Delta\theta_{2}\to
f_{r_{P_{12}}}}^{\\{\rho\\}}.$

Iterate - Outer Loop:
$\bm{\nu}\mathbf{=0,1,2,\dots};\bm{\rho}=\bm{\tau(\bm{\nu})}$

  1. 11.

Variable nodes compute marginals according to (18), e.g.:

$\begin{gathered}p(\Delta\theta_{2})\propto\mathcal{N}(\Delta\hat{\theta}_{2}|\Delta\theta_{2},\hat{\xi}_{\theta_{2}})\end{gathered}$
$\displaystyle\cfrac{1}{\hat{\xi}_{\Delta\theta_{2}}^{\\{\nu\\}}}$
$\displaystyle=\cfrac{1}{\xi_{\theta_{2}}}+\cfrac{1}{\xi_{f_{r_{P_{12}}}\to\Delta\theta_{2}}^{\\{\rho\\}}}+\cfrac{1}{\xi_{f_{r_{P_{23}}}\to\Delta\theta_{2}}^{\\{\rho\\}}}$
$\displaystyle\Delta\hat{\theta}_{2}^{\\{\nu\\}}$
$\displaystyle=\Bigg{(}\cfrac{r_{\theta_{2}}^{\\{\nu\\}}}{\xi_{\theta_{2}}}+\cfrac{r_{f_{r_{P_{12}}}\to\Delta\theta_{2}}^{\\{\rho\\}}}{\xi_{f_{r_{P_{12}}}\to\Delta\theta_{2}}^{\\{\rho\\}}}+\cfrac{r_{f_{r_{P_{23}}}\to\Delta\theta_{2}}^{\\{\rho\\}}}{\xi_{f_{r_{P_{23}}}\to\Delta\theta_{2}}^{\\{\rho\\}}}\Bigg{)}\hat{\xi}_{\Delta\theta_{2}}^{\\{\nu\\}}.$

  2. 12.

Variable nodes update the state variables, e.g.:

$\begin{gathered}\theta_{2}^{\\{\nu+1\\}}=\theta_{2}^{\\{\nu\\}}+\Delta\hat{\theta}_{2}^{\\{\nu\\}}.\end{gathered}$

  3. 13.

Repeat steps 4-13 until convergence111Note that, after step 8,, initialization
factor nodes are removed from the FG. Also in each iteration, virtual factor
nodes repeat the same message as in the initial step and messages from a
virtual factor node to a variable node should not be included in calculation
of the marginals..

##  V Numerical Results

The IEEE 14 bus test case shown in Fig. 6 is used to analyse performance of
the proposed algorithm. The set of measurements obtained from 61 measurement
devices that measure active and reactive power flow, active and reactive
injection power, bus voltage magnitude and bus voltage angle are selected so
that the system is observable.

Figure 6: The IEEE 14 bus test case

###  V-A Simulation setup

From a given IEEE 14 bus test case and the set of measurements, we define the
corresponding factor graph. Measurement values are generated using the AC
power flow analysis, additionally corrupted by Gaussian white noise of
variance $\sigma^{2}$. For each value of variance $\sigma^{2}$, using Monte
Carlo approach, we generate 1000 random sets of measurement values and feed
them to the proposed BP-based SE algorithm in order to obtain the average
performance results.

The convergence of the proposed algorithm is tracked by observing the root
mean square error (RMSE) after each outer iteration:

$\begin{gathered}\mathrm{RMSE}(\nu)={\cfrac{1}{n}||\mathbf{\hat{x}}-\mathbf{x}^{\nu}||_{2}},\end{gathered}$
(19)

where $\mathbf{\hat{x}}$ is the non-linear WLS solution, while
$\mathbf{x}^{\nu}$ represents the current iterate solution of the BP
algorithm.

We investigate two simulation scenarios using ”flat start” ($V_{i}=1$,
$\theta_{i}=0$, $i=1,\dots,N$). In the first scenario, we test the algorithm
convergence by measuring RMSE using the described Monte Carlo approach for
different values of measurement variances from the set $\sigma^{2}=$
$\\{\sigma_{1}^{2}$, $\sigma_{2}^{2}$, $\sigma_{3}^{2}$, $\sigma_{4}^{2}\\}$
$=\\{0.01^{2}$, $0.001^{2}$, $0.0001^{2}$, $0.00001^{2}\\}\;[\mbox{p.u.}]$.
For every outer iteration $\nu$, the number of inner iterations is defined as
$\tau(\nu)=\nu^{q}$, where $q$ is inner iteration exponent, which we set to
$q=4$. In the second scenario, we change the inner iteration number exponent
$q$ while keeping the variance fixed at the values $\sigma_{1}^{2}$ (high
noise level) and $\sigma_{4}^{2}$ (low noise level). The convergence is
discussed and presented in the following subsection.

###  V-B Simulation results

The set of measurements defines the topology of the FG and for almost all
placements of measurement devices of interest, the corresponding FG will have
loops222Note that, even if the physical power network has the radial structure
(tree structure), the FG will be loopy. An exception occurs, for example, for
the scenario of a radial network in which only end buses (and no internal
buses on a radial line) are allowed to contain power injection measurements..
It is well known that, in general, loopy BP does not converge to correct
marginals, e.g., specific inputs may lead to an oscillatory behaviour of
messages [12]. Based on extensive numerical studies, we presented in [13] a
heuristic solution to improve the convergence of the BP algorithm:

$\begin{gathered}\mu_{f\to x}^{\rho}=[1-\delta(p)]\cdot\mu_{f\to
x}^{\rho}+\delta(p)\cdot\alpha\cdot[\mu_{f\to x}^{\rho-1}+\mu_{f\to
x}^{\rho}],\end{gathered}$ (20)

where $\delta(p)\in\\{0,1\\}$ is a Bernoulli random variable with parameter
$p$, independently sampled for each message $\mu_{f\to x}^{\rho}$, and
$\alpha$ is the weighting coefficient. Namely, we modify updates of factor to
variable node messages in the inner iteration loop $\mu_{f\to x}^{\rho}$. For
the values of $p\in[0.4,0.6]$ and $\alpha=0.5$, our numerical studies showed
that the BP algorithm always converged successfully to the WLS solution.

Fig. 7 shows the convergence behaviour for the first simulation scenario
described in the previous subsection. The algorithm is terminated after
$\nu=7$ outer iterations with the total number of inner iterations
$\rho_{t}=\sum_{\nu=1}^{7}\nu^{4}=4676$. The figure demonstrates that the
proposed algorithm converges to WLS solution over a wide range of noise
variances. As expected, the convergence behaviour improves as the noise levels
decreases.

Figure 7: Convergence performance for different variances $\sigma^{2}$

The convergence behaviour for different values of $q$ for the set of
measurements with variances $\sigma_{1}^{2}$ (high noise level) and
$\sigma_{4}^{2}$ (low noise level) is shown in Fig. 8 and Fig. 9.

Figure 8: Convergence performance for different $q$ and $\sigma_{1}^{2}$

For different values of $q$, the number of outer iterations $\nu$ is selected
in such a way that the resulting number of inner iterations $\rho$ is
approximately the same. For both high and low noise levels, the value of inner
iteration number exponent $q=4$ is identified to provide fastest convergence
to the WLS solution for a given number of inner iterations. Note that for
insufficient value of exponent $q$ (e.g., $q=2$), the convergence speed of the
proposed algorithm will be dramatically reduced. The fastest convergence
behaviour is obtained for $q=4$ (although $q=5$ also shows good performance).
For too large $q$ value ($q\geq 6$), the convergence speed will decrease as
compared to $q=4$.

Figure 9: Convergence performance for different $q$ and $\sigma_{4}^{2}$

##  VI Conclusion

In this paper, we presented the BP solution of the AC SE problem that can be
interpreted as a fully distributed Gauss-Newton method. The proposed BP
algorithm converges to the same solution as the centralized WLS state
estimator. For the future work, we plan to compare the proposed algorithm in
the multi-area SE setup in terms of performance, convergence and computational
cost with the solutions recently proposed in the literature.

## Acknowledgment

This project has received funding from the EU 7th Framework Programme for
research, technological development and demonstration under grant agreement
no. 607774.

## APPENDIX
Measurement functions and Jacobian elements

The measurement functions
$\mathbf{h}(\mathbf{x})\equiv\mathbf{h}(\mathbf{V},\bm{\uptheta})$ that
connect measurements $\mathbf{z}$ to state variables
$\mathbf{x}\equiv(\bm{\uptheta},\mathbf{V})$ and corresponding Jacobian
elements are described below.

The measurement function for active power flow at the branch that connects
buses $i$ and $j$:

$\displaystyle h_{P_{ij}}(\cdot)$
$\displaystyle=V_{i}^{2}(g_{ij}+g_{si})-V_{i}V_{j}(g_{ij}\cos\theta_{ij}+b_{ij}\sin\theta_{ij}),$

where $V_{i}$ and $V_{j}$ are bus voltage magnitudes, while
$\theta_{ij}=\theta_{i}-\theta_{j}$ is the bus voltage angle difference
between buses $i$ and $j$. The parameters in above equations include the
conductance $g_{ij}$ and susceptance $b_{ij}$ of the branch, as well as the
conductance $g_{si}$ of the branch shunt element connected at the bus $i$. The
Jacobian expressions corresponding to $h_{P_{ij}}(\cdot)$ are as follows:

$\displaystyle\cfrac{\mathrm{\partial}{h_{P_{ij}}(\cdot)}}{\mathrm{\partial}\theta_{i}}$
$\displaystyle={V}_{i}{V}_{j}(g_{ij}\mbox{sin}\theta_{ij}-b_{ij}\mbox{cos}\theta_{ij})$
$\displaystyle\cfrac{\mathrm{\partial}{{h_{P_{ij}}}(\cdot)}}{\mathrm{\partial}\theta_{j}}$
$\displaystyle=-{V}_{i}{V}_{j}(g_{ij}\mbox{sin}\theta_{ij}-b_{ij}\mbox{cos}\theta_{ij})$
$\displaystyle\cfrac{\mathrm{\partial}{{h_{P_{ij}}}(\cdot)}}{\mathrm{\partial}V_{i}}$
$\displaystyle=-{V}_{j}(g_{ij}\mbox{cos}\theta_{ij}+b_{ij}\mbox{sin}\theta_{ij})+2V_{i}(g_{ij}+g_{si})$
$\displaystyle\cfrac{\mathrm{\partial}{h_{{P_{ij}}}(\cdot)}}{\mathrm{\partial}V_{j}}$
$\displaystyle=-{V}_{i}(g_{ij}\mbox{cos}\theta_{ij}+b_{ij}\mbox{sin}\theta_{ij}).$

The measurement function for reactive power flow at the branch that connects
buses $i$ and $j$:

$\displaystyle h_{Q_{ij}}(\cdot)$
$\displaystyle=-V_{i}^{2}(b_{ij}+b_{si})-V_{i}V_{j}(g_{ij}\sin\theta_{ij}-b_{ij}\cos\theta_{ij}),$

where $b_{si}$ is susceptance of the branch shunt element connected at the bus
$i$. The Jacobian expressions corresponding to $h_{Q_{ij}}(\cdot)$ are as
follows:

$\displaystyle\cfrac{\mathrm{\partial}{h_{{Q_{ij}}}(\cdot)}}{\mathrm{\partial}\theta_{i}}$
$\displaystyle=-{V}_{i}{V}_{j}(g_{ij}\mbox{cos}\theta_{ij}+b_{ij}\mbox{sin}\theta_{ij})$
$\displaystyle\cfrac{\mathrm{\partial}{h_{{Q_{ij}}}(\cdot)}}{\mathrm{\partial}\theta_{j}}$
$\displaystyle={V}_{i}{V}_{j}(g_{ij}\mbox{cos}\theta_{ij}+b_{ij}\mbox{sin}\theta_{ij})$
$\displaystyle\cfrac{\mathrm{\partial}{h_{{Q_{ij}}}(\cdot)}}{\mathrm{\partial}V_{i}}$
$\displaystyle=-{V}_{j}(g_{ij}\mbox{sin}\theta_{ij}-b_{ij}\mbox{cos}\theta_{ij})-2V_{i}(b_{ij}+b_{si})$
$\displaystyle\cfrac{\mathrm{\partial}{h_{{Q_{ij}}}(\cdot)}}{\mathrm{\partial}V_{j}}$
$\displaystyle=-{V}_{i}(g_{ij}\mbox{sin}\theta_{ij}-b_{ij}\mbox{cos}\theta_{ij}).$

The measurement function for active injection power into the bus $i$:

$\displaystyle h_{P_{i}}(\cdot)$
$\displaystyle=V_{i}\sum_{j\in\mathcal{H}}V_{j}(G_{ij}\cos\theta_{ij}+B_{ij}\sin\theta_{ij}),$

where $\mathcal{H}$ is the set of buses incident to the bus $i$, including the
bus $i$. The parameters $G_{ij}$ and $B_{ij}$ are conductance and susceptance
of the complex bus matrix. The Jacobian expressions corresponding to
$h_{P_{i}}(\cdot)$ are:

$\displaystyle\cfrac{\mathrm{\partial}{h_{P_{i}}(\cdot)}}{\mathrm{\partial}\theta_{i}}$
$\displaystyle={V}_{i}\sum_{j\in\mathcal{H}\setminus
i}{V}_{j}(-G_{ij}\mbox{sin}\theta_{ij}+B_{ij}\mbox{cos}\theta_{ij})$
$\displaystyle\cfrac{\mathrm{\partial}{h_{P_{i}}(\cdot)}}{\mathrm{\partial}\theta_{j}}$
$\displaystyle={V}_{i}{V}_{j}(G_{ij}\mbox{sin}\theta_{ij}-B_{ij}\mbox{cos}\theta_{ij})$
$\displaystyle\cfrac{\mathrm{\partial}{h_{P_{i}}(\cdot)}}{\mathrm{\partial}V_{i}}$
$\displaystyle=\sum_{j\in\mathcal{H}\setminus
i}{V}_{j}(G_{ij}\mbox{cos}\theta_{ij}+B_{ij}\mbox{sin}\theta_{ij})+2{V}_{i}G_{ii}$
$\displaystyle\cfrac{\mathrm{\partial}{h_{P_{i}}(\cdot)}}{\mathrm{\partial}V_{j}}$
$\displaystyle={V}_{i}(G_{ij}\mbox{cos}\theta_{ij}+B_{ij}\mbox{sin}\theta_{ij}),$

where $\mathcal{H}\setminus i$ is the set of buses incident to the bus $i$.

The measurement function for reactive injection power into the bus $i$:

$\displaystyle h_{Q_{i}}(\cdot)$
$\displaystyle=V_{i}\sum_{j\in\mathcal{H}}V_{j}(G_{ij}\sin\theta_{ij}-B_{ij}\cos\theta_{ij}),$

with Jacobian expressions:

$\displaystyle\cfrac{\mathrm{\partial}{h_{Q_{i}}(\cdot)}}{\mathrm{\partial}\theta_{i}}$
$\displaystyle={V}_{i}\sum_{j\in\mathcal{H}\setminus
i}{V}_{j}(G_{ij}\mbox{cos}\theta_{ij}+B_{ij}\mbox{sin}\theta_{ij})$
$\displaystyle\cfrac{\mathrm{\partial}{h_{Q_{i}}(\cdot)}}{\mathrm{\partial}\theta_{j}}$
$\displaystyle={V}_{i}{V}_{j}(-G_{ij}\mbox{cos}\theta_{ij}-B_{ij}\mbox{sin}\theta_{ij})$
$\displaystyle\cfrac{\mathrm{\partial}{h_{Q_{i}}(\cdot)}}{\mathrm{\partial}V_{i}}$
$\displaystyle=\sum_{j\in\mathcal{H}\setminus
i}{V}_{j}(G_{ij}\mbox{sni}\theta_{ij}-B_{ij}\mbox{cos}\theta_{ij})-2{V}_{i}B_{ii}$
$\displaystyle\cfrac{\mathrm{\partial}{h_{Q_{i}}(\cdot)}}{\mathrm{\partial}V_{j}}$
$\displaystyle={V}_{i}(G_{ij}\mbox{sin}\theta_{ij}-B_{ij}\mbox{cos}\theta_{ij}).$

The measurement function for current magnitude at the branch connecting buses
$i$ and $j$:

$\begin{gathered}h_{I_{ij}}(\cdot)=[aV_{i}^{2}+bV_{j}^{2}-2V_{i}V_{j}(c\cos\theta_{ij}-d\sin\theta_{ij})]^{1/2}\end{gathered}$
$\displaystyle a$ $\displaystyle=(g_{ij}+g_{si})^{2}+(b_{ij}+b_{si})^{2};$
$\displaystyle b$ $\displaystyle=g_{ij}^{2}+b_{ij}^{2}$ $\displaystyle c$
$\displaystyle=g_{ij}(g_{ij}+g_{si})+b_{ij}(b_{ij}+b_{si});$ $\displaystyle d$
$\displaystyle=g_{ij}b_{si}-b_{ij}g_{si}.$

The Jacobian expressions corresponding to current magnitude measurement
function $h_{I_{ij}}(\cdot)$ are as follows:

$\displaystyle\cfrac{\mathrm{\partial}{h_{I_{ij}}(\cdot)}}{\mathrm{\partial}\theta_{i}}$
$\displaystyle=\cfrac{V_{i}V_{j}(d\cos\theta_{ij}+c\sin\theta_{ij})}{I_{ij}}$
$\displaystyle\cfrac{\mathrm{\partial}{h_{I_{ij}}(\cdot)}}{\mathrm{\partial}\theta_{j}}$
$\displaystyle=-\cfrac{V_{i}V_{j}(d\cos\theta_{ij}+c\sin\theta_{ij})}{I_{ij}}$
$\displaystyle\cfrac{\mathrm{\partial}{h_{I_{ij}}(\cdot)}}{\mathrm{\partial}V_{i}}$
$\displaystyle=\cfrac{V_{j}(d\sin\theta_{ij}-c\cos\theta_{ij})+aV_{i}}{I_{ij}}$
$\displaystyle\cfrac{\mathrm{\partial}{h_{I_{ij}}(\cdot)}}{\mathrm{\partial}V_{j}}$
$\displaystyle=\cfrac{V_{i}(d\sin\theta_{ij}-c\cos\theta_{ij})+bV_{j}}{I_{ij}}.$

The Jacobian expressions corresponding to voltage magnitude and voltage angle
measurement functions are as follows:

$\displaystyle\cfrac{\mathrm{\partial}{{V_{i}}}}{\mathrm{\partial}V_{i}}=1;$
$\displaystyle\cfrac{\mathrm{\partial}{{\theta_{i}}}}{\mathrm{\partial}\theta_{i}}=1.$

## References

  * [1] A. Monticelli, _State Estimation in Electric Power Systems: A Generalized Approach_ , ser. Kluwer international series in engineering and computer science. Springer US, 1999.
  * [2] A. Abur and A. Expósito, _Power System State Estimation: Theory and Implementation_ , ser. Power Engineering. Taylor & Francis, 2004.
  * [3] F. F. Wu, K. Moslehi, and A. Bose, “Power system control centers: Past, present, and future,” _Proceedings of the IEEE_ , vol. 93, no. 11, pp. 1890–1908, Nov 2005.
  * [4] J. Pearl, _Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference_. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc., 1988.
  * [5] C. M. Bishop, _Pattern Recognition and Machine Learning_. Springer, 2006.
  * [6] Y. Weiss and W. T. Freeman, “On the optimality of solutions of the max-product belief-propagation algorithm in arbitrary graphs,” _IEEE Transactions on Information Theory_ , vol. 47, no. 2, pp. 736–744, Feb 2001.
  * [7] Y. Hu, A. Kuh, A. Kavcic, and D. Nakafuji, “Real-time state estimation on micro-grids,” in _Neural Networks (IJCNN), The 2011 International Joint Conference on_ , July 2011, pp. 1378–1385.
  * [8] Y. Weng, R. Negi, and M. Ilic, “Graphical model for state estimation in electric power systems,” in _Smart Grid Communications, 2013 IEEE International Conference on_ , Oct 2013, pp. 103–108.
  * [9] F. C. Schweppe and D. B. Rom, “Power system static-state estimation, part ii: Approximate model,” _IEEE Transactions on Power Apparatus and Systems_ , vol. PAS-89, no. 1, pp. 125–130, Jan 1970.
  * [10] P. Hansen, V. Pereyra, and G. Scherer, “Least squares data fitting with applications.” Johns Hopkins University Press, 2012, ch. 9, p. 166.
  * [11] H. A. Loeliger, J. Dauwels, J. Hu, S. Korl, L. Ping, and F. R. Kschischang, “The factor graph approach to model-based signal processing,” _Proc. of the IEEE_ , vol. 95, no. 6, pp. 1295–1322, 2007.
  * [12] K. P. Murphy, Y. Weiss, and M. I. Jordan, “Loopy belief propagation for approximate inference: An empirical study,” in _Proc. of the 15th Conference on Uncertainty in Artificial Intelligence_ , San Francisco, USA, 1999, pp. 467–475.
  * [13] M. Cosovic and D. Vukobratovic, “State estimation in electric power systems using belief propagation: An extended DC model,” in _Signal Processing Advances in Wireless Communications SPAWC, 2016. 17th IEEE Workshop_ , Edinburgh, United Kingdom, Jul. 2016.


# Bright Spatially Coherent Beam from Carbon Nanotube Fiber Field Emission
Cathode

Taha Y. Posos<EMAIL_ADDRESS>Electrical and Computer Engineering, Michigan
State University, MI 48824  Sergey V. Baryshev<EMAIL_ADDRESS>Electrical and
Computer Engineering, Michigan State University, MI 48824

###### Abstract

Large area carbon nanotube (CNT) cathodes made from yarns, films or fibers
have long been promising as next generation electron sources for high power
radio frequency (rf) and microwave vacuum electronic devices. However,
experimental evidence have been pointing out spatial incoherence of the
electron beam produced by such cathodes that, in turn, impeded the progress
toward high brightness CNT electron sources and their practical applications.
Indeed, typically large area CNT fibers, films or textiles emit stochastically
across their physical surface at large emission angles and with large
transverse spread, meaning large emittance and hence low brightness. In this
work, using high resolution field emission microscopy, we demonstrate that
conventional electroplating of hair-thick CNT fibers followed by a femtosecond
laser cutting, producing emitter surface, solves the described incoherent
emission issues extremely well. Strikingly, it was observed that the entire
(within the error margin) cathode surface of a radius of approximately
$75\,\mu\text{m}$ emitted uniformly (with no hot spots) in the direction of
the applied electric field. The normalized emittance on the fiber surface was
estimated of 52 nm with brightness of
$>$$10^{15}\frac{\text{A}}{\text{m}^{2}\text{rad}^{2}}$ (or $>$$10^{7}$ A
m-2sr-1V-1) estimated for pulsed mode operation.

field emission cathode, high beam brightness, low beam emittance

##  I Introduction

In early 1990’s, carbon nanotechnology revolution introduced a plethora of new
advanced materials among which the carbon nanotube was notoriously attractive
for making nanoscale field effects devices, including vacuum devices. Many
labs studied effects associated with field emission from a single CNT or
arrays with counted number of isolated CNTs.[1, 2, 3, 4, 5, 6] Control over
fabrication and emission of single CNT field emission devices was excellent
and many field emission devices were demonstrated, e.g. field emission
radio[7] or field emission transistor[8], amplifier[9], and many others.[5, 6]

Figure 1: Typical micrograph showing large beam transverse spread and
nonuniformity. The blue circle marks the cathode’s position behind the imaging
screen and its size.

In order to increase the output power, macrosocopic large area CNT fibers,
films, yarns and fabrics started to be used to increase the operating currents
from pico-/nano- to many amperes. Here, CNTs were thought to replace legacy
velvets.[10] The multiple benefits of CNT fibers over legacy technology are
low turn-on voltage and high emission current at relatively low operating
electric field due to inherent high field enhancement factor, and high
electrical and thermal conductivities.[11] It was conventionally assumed that
emission would be uniform, i.e. uniformity would translate from previously
studied arrays of counted CNTs to the large area CNT fibers. However, recent
studies that employed field emission microscopy illustrated that emission is
never uniform and moreover that the emission area is a function of the
electric field (making it cumbersome for calculating current densities.)
Fig.1, reproduced from our past work,[12] highlights other important issue of
the large transverse spread of the emitted beam where beam lands on the
imaging screen millimeters away from the physical location of the cathode
source (blue circle) after travelling only a millimeter between the cathode
and the anode. This clearly points out a very large emittance and therefore
very low brightness, making CNT fiber cathodes impractical for applications
like rf or microwave traveling wave tubes (operating in GHz range), microscopy
and bright X-ray sources for medicine or active scanning. Another issue arises
from that—because all the current emerges from a few active spots, it leads to
local heating, microbreakdowns[12] and short-lived cathodes.

After experimenting with many fiber arrangements, we found that it is tiny
singular fibrils (comprising braided fibers) that set loose due to thermal and
field related stress[12, 13] and that eventually focus the field due to their
high aspect ratio and become point-like randomized intense electron emitters
eventually exploding (seen as micro-breakdowns) and re-populating surrounding
areas with more new-born fibrils. This process repeats itself until the
cathodes stops operating while emission always look like a family of single
electron rays going in many directions that are not aligned with the desired
main longitudinal propagation direction, such as in Fig.1. To mitigate this
issue, hypothesised to be the major problem, we study new cathode production
technology where fibers are electroplated with Ni and laser cut; all to
suppress the fibril occurrence and regeneration. Through experimental
measurements and electrostatic and beam dynamics modeling, emission uniformity
and beam brightness were analyzed.

##  II Experimental

To prepare the field emission CNT fiber cathodes, a commercially available
CNTs fiber from DexMat, Inc was used. The fiber is made by a wet-spinning
technology [11]—pre-grown arrays of CNTs are dissolved in an acid to form a
spinnable liquid dope that is extruded through a spinneret into coagulant bath
to remove acid, and then dried in an oven. The resulting product is highly
aligned and densely packaged CNTs in a form of a fiber. DexMat fibers have
high electrical and thermal conductivity. Such fibers were shown to feature
anisotropic field emission[14] that is emission takes place along the fiber
(not from side walls) which is a great property allowing for control over
emittance. Raman spectroscopy shows the G peak positions at 1583 cm-1
suggesting rich crystalline graphitic content as expected from high quality
fibers (Fig.2A).

To mitigate the described stray fibril problem, few fibers of the described
kind were placed side by side and electroplated with Ni in an electrochemical
bath and flush cut from the top to the required length of about $5\,\text{mm}$
with a femtosecond micromachining laser beam. Then, it was welded on $1\times
1\,\text{inch}$ Ni base. The final fabricated structure can be seen in Fig.2B.

We tested two samples referred to as Sample A and Sample B through the rest of
the paper. In Sample A fibers were twisted and in Sample B were not, i.e.
simply place along each others. As scanning electron microscopy (SEM)
demonstrates in Fig.2C, additional fiber twisting enabled a dense core in
Sample A, while Sample B (Fig.2D) has visible voids between the individual
fibers. Otherwise, both sample have fiber core diameter of $\sim
150\,\mu\text{m}$ and Ni shell thickness of $\sim 50\,\mu\text{m}$. Sample A
and Sample B have height of $4.8\,\text{mm}$ and $4.6\,\text{mm}$ respectively
(see Fig.2B).

Figure 2: (A) Raman spectra of the cathode surface showing a crystalline
graphitic peak. (B) Electroplated CNT fiber welded on a Ni base;
$\text{H}=4.8\,\text{mm}$ for Sample A and $\text{H}=4.6\,\text{mm}$ for
Sample B. (C) and (D) SEM images of Sample A and Sample B, where scale bars
are 50 $\mu$m.

DC current tests and field emission microscopy were performed in our custom
field emission microscope described in great detail in Ref.15. Images were
processed by a custom image processing algorithm FEpic described
elsewhere.[16]

##  III Field Emission Imaging and Conditioning

After sample was installed and gap was tuned using a doublet of two orthogonal
optical microscopes, the physical location of the fiber is determined and
labeled. To do that, the test chamber is illuminated. Because the imaging
anode YAG screen is semitransparent, the location of the fiber can be
immediately seen and captured by photographing. The core of the fiber is
marked with red circle for reference in Figs.3A and 4A.

After that voltage is applied and field emission images are taken concurrently
with I-V curves. Fig.3B shows emission micrograph of the Sample A. The
improvement is immediately obvious when compared with Fig.1. First, the
emission spot appears exactly at the optical projection of the cathode. This
means beam divergence angle is small, so emittance can be expected to be low.
Second, there is only a single spot and its size is comparable to the size of
the fiber core—this is an indication of uniformity and small angular spread of
the electron beam. The same exact behaviour was observed for Sample B, as
given in Fig.4. No evidence suggesting the stray fibril issue was observed for
neither cathode. These results highlight that such a simple electroplating
strategy is extremely effective at yielding emission uniformity and spatial
coherence, thereby boosting the transverse beam brightness.

Figure 3: A) Sample A seen through the YAG screen when the lights is on in the
chamber. Its fiber core is marked with red circle. B) FE micrograph of the
same region at the gap of $200\,\mu\text{m}$. Figure 4: A) Sample B seen
through the YAG screen when the lights is on in the chamber. Its fiber core is
marked with red circle. B) FE micrograph of the same region at the gap of
$200\,\mu\text{m}$. Figure 5: A) The conditioning scheme of Sample A. Both
ramp up and down curves are shown. There is a clear decrease in performance.
B) Conditioning scheme of Sample B. Only ramp-up curves are shown. There is no
considerable change in performance.

Fig.5 shows cathode conditioning I-V curves for Sample A and B at
interelectrode gap of 200 $\mu$m. Conditioning[17, 18, 12] or cycling, where
the applied voltage is ramped up and down to a progressively higher number in
every consecutive cycle until the desired operating current is achieved, is a
crucial procedure to maximize field emission cathode performance and ensure
stability at the operating point.

We found that the electroplated fibers should be conditioned with small
incremental steps to avoid adverse effects such as sudden burn-down. Fig.5A
demonstrates the case where the maximal field was doubled with respect to the
previous conditioning cycle, such that the current went up from 3 to 22
$\mu$A. Next, after completing the ramp down Sample A stopped working
completely, which was possibly due to applying electric power that exceed that
in the previous run by more than an order of magnitude. Its operation could
not be rejuvenated by applying higher electric fields. This is unlike a
conditioning scheme that was used for Cathode B as shown in Fig.5B. The
emission current was doubled at every conditioning cycle: up to 1 $\mu$A and
down to 0, then to 2 $\mu$A, to 4 $\mu$A, to 8 $\mu$A, to 16 $\mu$A, and
finally to 32 $\mu$A. By doing so, Cathode B was conditioned $softly$
(compared to Sample A) maintaining and enhancing its performance: $i$) the
resulting operating field went up and doubled, reaching same exact value where
Cathode A burned down; $ii$) turn-on field and field enhancement factor
remained nearly the same meaning that Cathode B was conditioned to stably
sustain higher local field.

When it is compared to our past cathode designs, detailed in Ref.12, they emit
less at any given field. This is an expected result because (with stray
fibrils mitigated) the field enhancement is reduced. However, turn-on fields
are still very low, between 1 and 2.5 V/$\mu$m. Because the beam was tight
suggesting high current density we limited our measurements to between 10–100
$\mu$A as the power density deposition at the imaging screen could attain
above 1 kW/cm2 at the voltage source limit of 1100 V, thereby literally
drilling holes in it.[19, 20] At 1100 V, Sample A maxed out at 20 $\mu$A and
Sample B 30 $\mu$A, respectively. Again both cathodes had similar metrics.
Having these metrics and qualitative results in mind, a step was taken to
carry out more quantitative analysis and calculate cathodes’ emittance and
brightness. All detailed in the next section.

##  IV Emittance and Brightness

In the phase space $(x,x^{\prime})$, $x$ is spatial position and
$x^{\prime}=\frac{dx}{dz}=\frac{dx/dt}{dz/dt}=\frac{v_{x}}{v_{z}}$ is the
slope of the trajectory from longitudinal centrosymmetric axis of each
particle. Then, rms emittance $\tilde{\epsilon}_{x}$ is defined as

$\tilde{\epsilon}_{x}=\sqrt{\langle\Delta x^{2}\rangle\langle\Delta x^{\prime
2}\rangle-\langle\Delta x\Delta x^{\prime}\rangle}$ (1)

where $\Delta x=x-\langle x\rangle$ and $\Delta x^{\prime}=x^{\prime}-\langle
x^{\prime}\rangle$. For a beam with cylindrical symmetry in $(x,y)$ and
$(x^{\prime},y^{\prime})$ centered around zero, $\langle x\rangle$ and
$\langle x^{\prime}\rangle$ are zero. The cross-corelation term,
$\langle\Delta x\Delta x^{\prime}\rangle$ can be removed with proper beam
optics.[21] Then, Eq.1 reduces to

$\tilde{\epsilon}_{x}=\sqrt{\langle x^{2}\rangle\langle x^{\prime
2}\rangle}=\sigma_{x}\sqrt{\langle x^{\prime 2}\rangle}.$ (2)

Rms emittance is a function of the beam energy as $x^{\prime}$ is changing
under acceleration, and is not useful while comparing beam or beams at
different energies. On the other hand, from Liouville’s theorem, normalized
emittance is a conserved quantity under acceleration as long as the beam is
only subjected to conservative forces. The relation between rms emittance and
normalized emittance is given by

$\epsilon_{x}^{\text{N}}=\gamma\beta\tilde{\epsilon}_{x}$ (3)

where $\gamma=\frac{1}{\sqrt{1-\beta^{2}}}$ is the Lorentz factor and
$\beta=\frac{v}{\text{c}}\approx\frac{v_{z}}{\text{c}}$. In our case,
$\gamma\approx 1$ because energy is $<=$1 keV. Mean-transverse energy, MTE, is
$\frac{1}{2}\text{m}_{\text{e}}\langle v^{2}\rangle$, where
$v^{2}=v_{x}^{2}+v_{y}^{2}$. Because of the cylindrical symmetry in
$(v_{x},v_{y})$, $\text{MTE}\approx\frac{1}{2}\text{m}_{\text{e}}\langle
2v_{x}^{2}\rangle=\text{m}_{\text{e}}\langle v_{x}^{2}\rangle$. Then, after
substituting Eq.2, in terms of MTE, Eq.3 becomes

$\epsilon_{x}^{\text{N}}=\frac{v_{z}}{\text{c}}\cdot\sigma_{x}\cdot\sqrt{\langle\frac{v_{x}^{2}}{v_{z}^{2}}\rangle}=\sigma_{x}\cdot\sqrt{\frac{\text{MTE}}{\text{m}_{\text{e}}\text{c}^{2}}}$
(4)

Practically, the normalized emittance at the cathode surface is calculated as
follows. If the radius of the uniformly emitting surface of the cathode is
$r_{\text{i}}$, $\sigma_{x}\approx r_{\text{i}}$ can be taken. Moreover, at
the surface, MTE is due to statistical distribution of electrons inside the
cathode itself. So, it is intrinsic, and is further redefined as
$\text{MTE}_{\text{i}}$. Then, Eq.4 becomes

$\epsilon_{x}^{\text{N}}=r_{\text{i}}\cdot\sqrt{\frac{\text{MTE}_{\text{i}}}{\text{m}_{\text{e}}\text{c}^{2}}}$
(5)

By using normalized emittance, normalized transverse brightness,
$B_{\text{N}}$, can be calculated as

$B_{\text{N}}=\frac{2\,I}{\epsilon_{x}^{\text{N}}\epsilon_{y}^{\text{N}}}$ (6)

where $I$ is the emitted current.
$\epsilon_{x}^{\text{N}}=\epsilon_{y}^{\text{N}}$ can be taken in cylindrical
symmetry.

Figure 6: Series of field emission micrographs of Sample B with the screen
moved away progressively from $200\,\mu\text{m}$ gap to $1600\,\mu\text{m}$
gap with $200\,\mu\text{m}$ steps. At each step, a micrograph is captured.

To determine $r_{\text{i}}$ and $\text{MTE}_{\text{i}}$ in Eq.5, we conducted
experimental measurements in combination with beam dynamics in GPT (General
Particle Tracer).[22] In the measurements, the imaging screen is moved away
from the cathode progressively. The voltage is set accordingly to keep the
current constant at 20 nA to enable strong beam image signal but avoid
additional beam expansion due to vacuum space charge effect. A micrograph at
each step is recorded (see Fig.6). Increase in the spot size due to larger
time of flight is measured. As it is seen in Fig.6, the spots are Gaussian in
nature with dense center and faint tails. Each spot can be modeled
mathematically with a cylindrically symmetric Gaussian as

$p=A\cdot\text{exp}\left(-\frac{(x-x_{c})^{2}+(y-y_{c})^{2}}{2\sigma_{\text{spot}}^{2}}\right)+C$
(7)

to extract projected transverse beam size. Here, $A$ is the amplitude,
$\sigma_{\text{spot}}$ is the standard deviation, $C$ is the background
offset, $p$ is the intensity, $(x,y)$ are the space dimensions,
$(x_{c},y_{c})$ are the coordinates of the peak.[16] The model parameters $A$,
$\sigma_{\text{spot}}$, and $C$ for each spot are computed with least-square
fitting method. After fitting, the emission spot diameter is taken as
$2\sigma_{\text{spot}}$. An exemplary 3D fitting done by FEpic for the beam
imaged at 200 $\mu$m gap is presented in Fig.7; with the black mesh surface
being the fitting surface. The resulting dependence of the spot size, first
measured (Fig.6) and processed by FEpic, versus distance is shown in Fig.8
with black solid circles. The data in the figure is only for Sample B. Because
Sample A burned down, studies similar to those presented in Fig.6 could not be
carried out.

Figure 7: The color surface shows the beam spot for $200\,\mu\text{m}$ gap in
3D. The black mesh surface shows its mathematical fit in Eq.7. Figure 8:
Comparison of experimental and computational final beam spot size as the
screen is moving away from the cathode. $\text{MTE}_{\text{i}}$ of
$250\,\text{eV}$ and $r_{\text{i}}$ of $75\,\mu\text{m}$ were used in GPT
modeling. Figure 9: Illustration of electric field computed in COMSOL for a
$1\,\text{mm}$ gap. The color plot shows the field magnitude and contour. The
dark region is the fiber core, and the gray region is the Ni shell. Figure 10:
In GPT: (A) Initial uniform beam distribution at the cathode surface in real-
space, where $r_{\text{i}}$ is the radius of the beam. (B) The initial
distribution in momentum-space, where $\beta_{x}=v_{x}/c$,
$\beta_{y}=v_{y}/c$, and $r_{\beta}$ is the radius. Final distribution in (C)
real space and (D) momentum space when the screen at is $1\,\text{mm}$.

To calculate the phase space volume for emittance and brightness estimations
and obtain $r_{\text{i}}$ and $\text{MTE}_{\text{i}}$, we switched to beam
dynamics in GPT by comparing it with the experiment in Fig.6. To do that, a
field map for each interelectrode gap was computed in COMSOL by solving
Poisson’s equation with given boundary conditions. An exemplary field
distribution for 1 mm gap is shown in Fig.9. Then, the field maps were
imported to GPT. In GPT, the initial particle distribution in the real space
$(x,y)$ (Fig.10A) and momentum space
$(\beta_{x},\beta_{y})=(v_{x}/\text{c},v_{y}/\text{c})$ (Fig.10B) has to be
set to propagate the beam. Because the cathode itself circular, we used a
circular uniform distribution with the radius $r_{\text{i}}$ in position-space
and the radius $r_{\beta}$ in $\beta$-space. $\text{MTE}_{\text{i}}$ of the
fiber is expected to be $250\,\text{meV}$.[23] To be used in GPT,
$\text{MTE}_{\text{i}}$ is converted into $r_{\beta}$ as

$r_{\beta}=\sqrt{\frac{4\,\text{MTE}_{\text{i}}}{\text{m}_{\text{e}}\text{c}^{2}}}.$
(8)

$\text{MTE}_{\text{i}}$ of 250 meV translates to $r_{\beta}$ of $1.4\times
10^{-3}$. At the same time, $r_{\text{i}}$ remains a free model parameter to
be found by finding the best agreement between GPT with the experiment.

The beam was launched at a charge corresponding to 20 nA and allowed to drift
through the distance corresponding to a specified interelectrode gap. Here,
the beam dynamics is computed self-consistently taking the COMSOL calculated
field. The final distribution in $(x,y)$ and $(\beta_{x},\beta_{y})$ was
captured at a distance corresponding to the imaging YAG screen of the
microscope and are shown in Fig.10C and D. In the $(x,y)$ space, the standard
deviation $\sigma_{gpt}$ of such projections were calculated for every
cathode-anode gap, and $2\sigma_{gpt}$ was taken as the computed beam diameter
(analogous to FEpic image processing of the experimental images). It was
established that in GPT, when $r_{\text{i}}$ was set to the physical radius of
the fiber core of $75\,\mu\text{m}$, and $r_{\beta}$ was set by MTE of 250
meV,[23] then the final diameter (spot size $2\sigma_{gpt}$) of the resulting
beam projection was in a very good quantitative agreement with the experiments
(see Fig.8). Note, the GPT results were fairly insensitive to MTE values set
between 25 (typical Fermi level value for CNT) and 250 meV, and magnification
was due to the radial field distribution. This points out that the divergence
between the experiment and GPT (seen for the gap ranging between 1 and 1.6 mm)
stems from the difference between the idealized computed and the actual field
distribution in the gap. The summary of the results in Fig.8 confirms that the
entire whole fiber surface actively and uniformly emitting with a small MTE.

Finally, using Eq.5 and substituting $r_{\text{i}}=75\,\mu\text{m}$ and
$\text{MTE}_{\text{i}}=250\,\text{eV}$, the upper limit of the normalized
emittance on the fiber cathode surface can be estimated as

$\epsilon_{x}^{N}=0.052\,\text{mm}\,\text{mrad}=52\,\text{nm}.$ (9)

From this, taking the measured current (limited to 10–100 $\mu$A due to
extremely high power density), as shown in Fig.5, the normalized brightness
for 50 $\mu$A dc current is $B_{N}=3.7\times
10^{10}\frac{\text{A}}{\text{m}^{2}\text{rad}^{2}}$. The same very fibers can
draw currents of 1–10 A when operated in pulsed mode with a pulse length of
100–300 ns.[24] Using the estimated emittance of $52\,\text{nm\,rad}$, the
brightness in the pulsed mode, preferable mode in most VED HPM applications,
attains a notable value of $B_{N}=4.4\times
10^{15}\frac{\text{A}}{\text{m}^{2}\text{rad}^{2}}$. This number is
outstanding and is comparable with brightness metrics in the the state-of-the-
art microwave/rf accelerator injectors.[25] This (phase space) brightness can
be converted into geometrical brightness, a definition of brightness commonly
employed in the electron microscopy literature. The geometrical reduced
brightness is defined as
$B_{r}^{G}=\frac{\text{d}I}{\text{d}\Omega}\frac{1}{U}\frac{1}{S_{\text{cathode}}}$,
where $\Omega$ is the solid angle, $U$ is the voltage at which the current $I$
is measured, and $S_{\text{cathode}}$ is the emission area of the cathode. Our
calculations show that in pulsed mode it could attain $B_{r}^{G}=5.7\times
10^{7}$ A m-2sr-1V-1. This number is within the range obtained for single CNT
emitters.[26]

##  V Conclusion

In conclusion, we presented a simple and efficient field emission cathode
design where CNT fiber core was plated with nickel shell. This design had two
important functions. First, it compresses the core, provides mechanical
strength thereby preventing stray fibril formation during conditioning and
operation. Second, such design (while slightly reducing field enhancement and
increasing turn-on field) reduces the fringing field on the CNT fiber and
therefore the defocusing radial field.

As field emission microscopy directly demonstrated, both tested cathodes
featured excellent spatially coherent emission. Field emission microscopy
aided by image processing and beam dynamics simulations confirmed that the
entire fiber core of 150 $\mu$m in diameter actively and uniformly emitted
electrons, as well as enabled phase space analysis. All of these combined
allowed to quantify the observed emission coherence through calculating
emittance and brightness. The extremely low emittance resulting in record
brightness highlight a simple and practical path forward for the CNT fiber
technology that has long been expected to advance high frequency vacuum power
devices but had limited success due to low brightness.

Finally and most importantly, it was demonstrated that the nanoscopic single
CNT cathode technology can be translated to the macroscopic fiber CNT level in
terms of emission uniformity. In other words, spatial coherence and uniformity
(intrinsic to a single CNT emitter) can be achieved in a CNT fiber comprised
out of billions of single CNT’s. The obtained brightness figures of merit
further confirm this technology translation in that ultimate single CNT
emitter brightness is feasible to attain for CNT fiber cathodes.

## References

  * Dean and Chalamala [1999] K. A. Dean and B. R. Chalamala, Journal of Applied Physics 85, 3832 (1999).
  * 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).
  * De Jonge  _et al._ [2002] N. De Jonge, Y. Lamy, K. Schoots, and T. H. Oosterkamp, Nature 420, 393 (2002).
  * Purcell  _et al._ [2002] S. T. Purcell, P. Vincent, C. Journet, and V. T. Binh, Phys. Rev. Lett. 88, 105502 (2002).
  * Teo  _et al._ [2004] K. Teo, R. Lacerda, M. Yang, A. Teh, L. Robinson, S. Dalal, N. Rupesinghe, M. Chhowalla, S. Lee, D. Jefferson, _et al._ , IEE Proceedings-Circuits, Devices and Systems 151, 443 (2004).
  * Li  _et al._ [2015] Y. Li, Y. Sun, and J. Yeow, Nanotechnology 26, 242001 (2015).
  * Jensen  _et al._ [2007] K. Jensen, J. Weldon, H. Garcia, and A. Zettl, Nano letters 7, 3508 (2007).
  * Wong  _et al._ [2005] Y. Wong, W. Kang, J. Davidson, B. Choi, W. Hofmeister, and J. Huang, Diamond and related materials 14, 2069 (2005).
  * Wong  _et al._ [2009] Y. M. Wong, W. P. Kang, J. L. Davidson, D. V. Kerns, J. Huang, and K. F. Galloway, IEEE transactions on electron devices 56, 738 (2009).
  * Shiffler  _et al._ [2012] D. Shiffler, S. Fairchild, W. Tang, B. Maruyama, K. Golby, M. LaCour, M. Pasquali, and N. Lockwood, IEEE Transactions on Plasma Science 40, 1871 (2012).
  * Behabtu  _et al._ [2013] N. Behabtu, C. C. Young, D. E. Tsentalovich, O. Kleinerman, X. Wang, A. W. K. Ma, E. A. Bengio, R. F. ter Waarbeek, J. J. de Jong, R. E. Hoogerwerf, S. B. Fairchild, J. B. Ferguson, B. Maruyama, J. Kono, Y. Talmon, Y. Cohen, M. J. Otto, and M. Pasquali, Science 339, 182 (2013).
  * 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).
  * Fairchild  _et al._ [2015] S. B. Fairchild, J. Boeckl, T. C. Back, J. B. Ferguson, H. Koerner, P. T. Murray, B. Maruyama, M. A. Lange, M. M. Cahay, N. Behabtu, C. C. Young, M. Pasquali, N. P. Lockwood, K. L. Averett, G. Gruen, and D. E. Tsentalovich, Nanotechnology 26, 105706 (2015).
  * Fairchild  _et al._ [2021] S. B. Fairchild, T. A. de Assis, J. H. Park, M. Cahay, J. Bulmer, D. E. Tsentalovich, Y. S. Ang, L. K. Ang, J. Ludwick, T. C. Back, and M. Pasquali, Journal of Applied Physics 129, 125103 (2021).
  * Baturin and Baryshev [2017] S. S. Baturin and S. V. Baryshev, Review of Scientific Instruments 88, 033701 (2017).
  * Posos  _et al._ [2021] T. Y. Posos, O. Chubenko, and S. V. Baryshev, ACS Applied Electronic Materials 3, 4990 (2021).
  * Millikan and Eyring [1926] R. A. Millikan and C. F. Eyring, Phys. Rev. 27, 51 (1926).
  * 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).
  * Jasim  _et al._ [1992] K. M. Jasim, R. Rawlings, and D. West, Journal of materials science 27, 1937 (1992).
  * Klimov  _et al._ [2016] A. Klimov, A. Zenin, Y. Yushkov, and E. Oks, in  _AIP Conference Proceedings_ , Vol. 1772 (2016) p. 040003.
  * Jarvis [2009] J. D. Jarvis,  _Development of high-brightness electron sources for free-electron lasers_ , Ph.D. thesis, Vanderbilt University (2009).
  * van der Geer and de Loos [2022] S. van der Geer and M. de Loos, General particle tracer (pulsar physics) (2022).
  * Back  _et al._ [2019] T. C. Back, G. Gruen, J. Park, P. T. Murray, J. Ludwick, M. Cahay, and S. B. Fairchild, AIP Advances 9, 065319 (2019).
  * Fairchild  _et al._ [2020] S. B. Fairchild, P. T. Murray, S. Portillo, and G. Dion, in  _2020 IEEE 21st International Conference on Vacuum Electronics (IVEC)_ (2020) pp. 19–20.
  * Rosenzweig  _et al._ [2019] J. B. Rosenzweig, A. Cahill, V. Dolgashev, C. Emma, A. Fukasawa, R. Li, C. Limborg, J. Maxson, P. Musumeci, A. Nause, R. Pakter, R. Pompili, R. Roussel, B. Spataro, and S. Tantawi, Phys. Rev. Accel. Beams 22, 023403 (2019).
  * de Jonge  _et al._ [2005] N. de Jonge, M. Allioux, J. T. Oostveen, K. B. K. Teo, and W. I. Milne, Phys. Rev. Lett. 94, 186807 (2005).